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Abstract 


Several  new  quadrature  sets  for  use  in  the  discrete  ordinates  method  of 
solving  the  Boltzmann  neutral  particle  transport  equation  are  derived.  These 
symmetric  quadratures  extend  the  traditional  symmetric  quadratures  by 
allowing  ordinates  perpendicular  to  one  or  two  of  the  coordinate  axes. 
Comparable  accuracy  with  fewer  required  ordinates  is  obtained. 

Quadratures  up  to  seventh  order  are  presented.  The  validity  and  efficiency  of 
the  quadratures  is  then  tested  and  compared  with  the  LQn  level  symmetric 
quadratures  relative  to  a  Monte  Carlo  benchmark  solution.  The  criteria  for 
comparison  include  current  through  the  surface,  scalar  flux  at  the  surface, 
volume  average  scalar  flux,  and  time  required  for  convergence.  Appreciable 
computational  cost  was  saved  when  used  in  an  unstructured  tetrahedral  cell 
code  using  highly  accurate  characteristic  methods.  However,  no  appreciable 
savings  in  computation  time  was  found  using  the  new  quadratures  compared 
with  traditional  Sn  methods  on  a  regular  Cartesian  mesh  using  the  standard 
diamond  difference  method.  These  quadratures  are  recommended  for  use  in 
three-dimensional  calculations  on  an  unstructured  mesh. 
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I.  Introduction 


This  research  developed  a  set  of  angular  quadratures  that  are 
computationally  efficient  when  used  with  the  discrete  ordinates  method  to 
solve  the  three-dimensional  Boltzmann  neutral  particle  transport  equation. 
The  quadrature  sets  contain  directions  perpendicular  to  one  or  more  cardinal 
directions.  These  sets  are  tested  for  accuracy  and  computational  efficiency. 
Performance  comparisons  are  made  with  traditional  (level  symmetric) 
quadrature  sets.  Types  of  problems  and  conditions  where  these  quadratures 
are  most  applicable  are  discussed. 

Background 

The  foundation  of  transport  theory  is  the  Boltzmann  transport 
equation  (BTE).  This  equation,  formulated  over  a  century  ago,  was  originally 
developed  for  the  study  of  the  kinetic  theory  of  gases  (1:  1).  Preliminary 
study  in  the  field  was  primarily  of  the  diffusion  of  light  by  the  atmosphere. 
Then  study  began  in  the  early  part  of  the  twentieth  century  on  investigating 
the  diffusion  of  energy  through  the  atmosphere  of  a  star  (2:  1).  The  scale  of 
these  problems  are  such  that  they  can  be  modeled  as  semi-infinite  media 
with  one-dimensional  geometry  and  therefore  the  methods  of  their  solution 
are  of  limited  application  (1:  1). 
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In  the  1940s,  interest  in  the  military  and  industrial  application  of 
nuclear  energy  stimulated  a  tremendous  amount  of  research  into  neutral 
particle  transport.  The  incredible  urgency  involved  in  the  research  of  nuclear 
energy  induced  by  the  events  of  World  War  II  necessarily  resulted  in  the 
development  and  use  of  approximate  methods  for  solving  the  linearized 
transport  equation  (2:1). 

Motivation 

Military  and  industrial  research  into  the  use  of  nuclear  energy  using 
actual  nuclear  material  has  decreased  substantially  in  recent  years.  The 
comprehensive  nuclear  test  ban  treaty  (CTBT),  if  ratified,  will  eliminate  our 
ability  to  obtain  any  further  real  data  on  new  weapons  designs  or  systems 
survivability  in  or  near  a  real  threat  environment.  Some  aspects  of  the 
radiation  environment  resulting  from  a  nuclear  detonation  are  simulated  at 
various  test  sites.  These  tests  can  only  approximate  the  actual  post 
detonation  environment  and  are  very  costly  (17).  It  is  currently  politically 
undesirable  for  private  industry  to  perform  research  using  significant  of 
amounts  nuclear  material  or  to  build  new  nuclear  research  facilities.  These 
and  other  factors  have  greatly  increased  the  need  for  accurate  computer 
modeling  of  nuclear  material  and  effects.  The  high  performance  computers 
needed  for  this  modeling  are  very  expensive,  and  therefore  anything  that  can 
increase  the  efficiency  of  these  machines  will  translate  directly  into 
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substantial  savings  by  increasing  their  productivity  and  extending  their 
useful  life.  Reduction  in  computation  time  also  increases  the  practicality  of 
analyzing  a  large  variety  of  similar  scenarios  for  threat  analysis  or  system 
optimization. 

The  Boltzmann  Transport  Equation 

The  Boltzmann  transport  equation  (Equation  1-1)  is  a  conservation 
equation  for  the  flux  of  neutral  particles  (1:  24).  The  particles  can  be 
neutrons,  photons,  or  any  other  neutral  particle  given  the  nuclear  data.  The 
angular  flux,  v| / ,  is  dependent  on  position  (r),  direction  of  motion  (Q),  on 
speed  or  energy  (v,  E),  and  time  (t).  Equation  (1-1)  represents  a  balance 
between  the  loss  rate  (right  side)  and  gain  rate  (left  side)  of  particles  that 
exist  at  each  point  of  this  seven  dimensional  phase  space  (3: 1-2); 

-  +  Q  •  V  +  at  (?,  E>  t)Ur,  E,  Q,  t)  = 

|_v  dt  J  (1-1) 

jdE  Jd0as  (r,  E  -►  Ef,  Q  •  Of,  t)q/(f,  V,  Q,  t) + s(f ,  E,  Q,  t) 

where  the  variables  are  defined  in  Table  1. 
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Table  1-1:  Variable  Definitions  for  the  BTE,  Equation  (1-1) 


Variable 

Description 

V 

magnitude  of  the  velocity 

t 

time 

r 

position 

E 

energy 

total  macroscopic  cross-section  for 
interaction  (  absorption  and  scatter  ) 

°s 

macroscopic  scattering  cross-section 

s 

total  non-scattering  source 

V 

angular  flux:  a  distribution  function 
of  particles  at  point  r,  with  energy  E, 
moving  in  direction  Q  at  time  t 

Q 

unit  vector  aligned  along  the 
streaming  direction  of  particles:  it  is 
often  shown  as  three  components,  or 
direction  cosines,  p,  p,  £,  defined  by 

P  —  Q  *  5  h  ^  ^  ^  *  ^2 

as  shown  in  Figure  1-1 

Figure  1-1:  Direction  Cosines 


Though  discrete  ordinates  is  valid  for  time  dependent  problems,  this 
treatment  will  assume  steady  state  conditions  with  all  time  dependence 
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suppressed;  therefore  the  first  term  in  Equation  (1-1)  (representing  the  time 
rate  of  change  of  particles  in  the  phase  space)  vanishes.  The  remaining  terms 
represent  the  steady  state  balance  equation.  The  second  term  in  the  brackets 
on  the  left  of  Equation  (1-1)  is  the  streaming  operator  and  represents  the  loss 
rate  due  to  particle  divergence.  The  final  term  in  the  brackets  on  the  left  is 
the  collision  operator  and  represents  the  loss  rate  due  to  particle  interaction 
with  the  medium.  This  interaction  could  be  absorption  (destroying  the 
particle),  or  scatter  (changing  the  particle’s  energy  or  direction)  (3: 1-2). 

The  first  term  on  the  right  of  Equation  (1-1)  is  the  gain  rate  due  to 
particles  traveling  in  other  directions  that  scatter  into  the  given  direction,  Q . 
As  a  consequence  of  the  isotropic  material  assumption,  the  distribution  of 
particles  scattering  into  a  given  direction  is  not  a  function  of  the  incident 
angle,  Q' ;  however  it  is  a  function  of  the  angle  between  the  incident  direction 
and  final  direction,  (q  •  Q').  Despite  the  isotropic  material  assumption,  this 
dependence  on  scattering  angle  means  that  scattering  may  be  anisotropic  (1). 
The  final  term  represents  the  gain  rate  from  production  of  particles  by  any 
source  mechanism.  The  source  can  be  internal,  such  as  radioactive  decay  and 
fission  or  external  such  as  solar  x-rays  or  an  incident  beam.  A  more  detailed 
discussion  of  the  BTE  and  definition  of  its  components  is  presented  in  section 
II. 

The  macroscopic  cross  sections  are  functions  of  position,  through 
spatially  varying  number  density  or  changes  in  material,  and  of  energy 
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through  the  microscopic  cross  sections  in  which  fundamental  properties  of  the 
isotopes  are  involved.  In  general  the  microscopic  cross  sections  could  have 
some  additional  implicit  spatial  variation.  For  example,  this  could  be  due  to 
a  temperature-dependence  through  Doppler  broadening  (1:  6).  This  research 
assumes  each  material  is  uniform  and  isotropic  and  therefore  no  such  implicit 
spatial  dependence  exists. 

The  physical  quantities  most  often  of  interest  to  be  found  via  the  BTE 
are  the  scalar  flux  <j)  (zeroth  angular  moment  of  the  vector  flux)  and  the 
vector  current  J  (first  angular  moment).  The  scalar  flux  is  of  interest  as  it 
represents  the  total  expected  particle  path  length  traveled  per  unit  volume  at 
a  given  location  in  the  medium.  This  will  determine  the  reaction  rates  for 
such  things  as  fission  and  neutron  activation.  The  vector  current  will 
determine  the  leakage  rate  from  one  region  to  the  next  or  through  boundaries 
(3: 1-3). 

Only  a  limited  number  of  analytic  and  semi-analytic  solutions  exist  for 
the  BTE.  Most  of  these  solutions  are  for  highly  idealized  problems.  Many 
ingenious  methods  such  as  discrete  ordinates,  Monte  Carlo,  even-parity, 
finite-elements,  and  Green  functions  have  been  developed  to  solve  the 
transport  equation  and  to  extend  the  application  of  such  knowledge.  The 
diffusion  equation  can  be  derived  from  the  BTE  using  several  simplifying 
assumptions  regarding  the  angular  dependence.  For  problems  with  nearly 
isotropic  scatter,  this  can  yield  approximate  results  (16:  2).  Of  particular 
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interest  are  methods  capable  of  providing  solutions  to  the  broad  range  of 
geometrical  configurations  found  in  nuclear  reactor  and  radiation  shielding 
applications  (1:  1).  From  a  military  perspective,  in  addition  to  the  application 
to  power  generation,  accurate  modeling  of  neutron  flux  allows  for  more 
precise  modeling  of  the  yield  and  effects  of  nuclear  weapons. 

The  advent  of  high-speed  computing  and  large  storage  capacity  has  led 
to  the  refinement  and  use  of  two  primary  numerical  methods  of  attaining  a 
solution  to  the  BTE:  the  method  of  discrete  ordinates  and  Monte  Carlo  (1:  2). 
Other  less  used  methods  will  not  be  examined.  The  theory  and  development 
of  the  Monte  Carlo  method  will  not  be  discussed  in  detail;  however,  Monte 
Carlo  solutions  will  be  used  as  benchmarks. 

Since  its  evolution  from  the  angular  segmentation  formulation  by 
Carlson  in  1957  discrete  ordinates  has  become  a  widely  used  method  for 
solving  the  integrodifferential  form  of  the  transport  equation.  The  discrete 
ordinates  method  involves  enforcing  the  transport  equation  only  at  discrete 
angular  directions  called  ordinates.  These  ordinates  are  selected  such  that 
the  flux  moments  may  be  evaluated  accurately  by  a  weighted  sum  (1:  118). 
For  example  the  scalar  flux  (zeroth  flux  moment)  is 

<j)(f,E,t)=  JdQv|/(r,E,Qt)«  ^wn\|/(r,E,Qn,t). 
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A  more  detailed  discussion  of  the  flux  moments  is  found  in  Chapter  II.  The 
advantages  of  the  discrete  ordinates  method  are  the  relatively  simple 
derivation  and  the  subsequent  ease  of  transformation  into  algorithms  of  good 
computational  efficiency  (1:  116).  It  also  lends  itself  well  to  the  discretizing  of 
energy  into  multiple  energy  groups.  A  distinct  advantage  over  the  Monte 
Carlo  method  is  that  it  provides  flux  and  current  data  everywhere  in  the 
problem  rather  than  only  at  a  limited  number  of  locations.  A  quadrature  set 
is  the  combination  of  discrete  angles  and  weights  used  in  a  weighted  sum  to 
evaluate  the  flux  moments. 

There  are  two  independent  angular  directions  for  Q .  The  directions 
are  parameterized  by  three  direction  cosines  that  obey  the  relationship 

p2  +  r|2  +  ^2  =  1  (1-3) 

where  p,  r\,  and  %  are  shown  in  Figure  1-1. 

Once  the  angular  approximation  has  been  made,  a  spatial 
discretization  scheme  must  be  used.  Computational  cost  and  storage 
requirements  are  directly  proportional  to  the  number  of  spatial  cells  and 
discrete  ordinates  used.  A  large  number  of  spatial  schemes  have  been 
formulated  for  use  in  discrete  ordinates  calculations.  They  include  linear 
methods  such  as  diamond  difference,  linear  discontinuous,  and  linear 
characteristic  and  non-linear  methods  such  exponential  characteristic  (9). 
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The  most  serious  drawback  to  the  discrete  ordinates  method  is  the 
buildup  of  truncation  errors  due  to  the  discretized  angular  and  spatial 
representations  (1:  131).  The  truncation  errors  can  result  in  random  error, 
which  limits  the  accuracy  of  the  results,  or  may  lead  to  physically  unrealistic 
results  such  as  negative  fluxes  or  sources.  Systematic  truncation  error  may 
lead  to  what  are  known  as  ray  effects  (21,  8,  3).  These  are  errors  caused  by 
the  discrete  ordinates  method  of  limiting  particle  motion  to  discrete 
directions  or  rays.  Flux  due  to  unscattered  particles  will  only  be  found  to 
occur  at  points  where  a  line  can  be  drawn  from  a  source  to  the  point  in  the 
direction  of  a  discrete  ordinate.  This  causes  the  scalar  flux  to  be  calculated 
higher  than  expected  at  points  along  discrete  ordinate  directions  and  lower 
between.  The  method  is  not  well  suited  to  geometry  with  a  strongly  peaked 
flux  in  a  given  direction.  Generally,  a  separate  method  must  be  used  to 
calculate  the  first  scatter  source  for  such  a  problem.  In  order  to  increase 
accuracy  of  the  discrete  ordinates  method  and  minimize  the  negative 
consequences,  it  is  either  necessary  to  increase  the  number  of  directions  in 
the  angular  quadrature,  thus  increasing  the  computation  time  and  storage 
requirements  of  the  computer  system,  or  to  develop  an  alternative 
quadrature  that  produces  less  error  with  fewer  directions.  This  research 
concentrates  on  increasing  the  accuracy  of  the  discrete  ordinates 
approximation  while  minimizing  the  number  of  angles. 
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The  primary  drawback  of  allowing  motion  perpendicular  to  one  or 
more  cardinal  directions  in  a  quadrature  set  is  that  mathematical  instability 
may  result  when  using  current  computer  codes.  The  resulting  zero 
components  of  flux  often  generate  run-time  errors.  The  advantages  of  using 
directions  parallel  to  cell  boundaries  are  (due  to  one  or  two  of  the  direction 
cosines  being  zero)  the  discrete  ordinates  equations  simplify  significantly  and 
fewer  directions  are  required  for  the  same  order  of  anisotropy,  thus  allowing 
for  increased  computational  efficiency. 

Except  for  the  simple  case  of  isotropic  scatter,  the  cross  section  for 
scattering  will  be  a  function  of  the  scattering  angle  as  well  as  energy . 
Separation  of  the  angular  and  energy  dependence  is  assumed.  Traditionally, 
cross  sections  are  then  expanded  in  orthogonal  Legendre  polynomials  (1:  13). 
The  order  of  this  expansion  is  another  limit  to  the  accuracy  obtainable  by  the 
discrete  ordinates  method. 

The  high  performance  computers  needed  to  perform  these  calculations 
are  very  expensive,  and  therefore  anything  that  can  increase  the  efficiency  of 
these  machines  will  translate  directly  into  substantial  savings  by  increasing 
their  useful  life.  Reduction  of  computation  time  also  increases  the 
practicality  of  analyzing  a  large  variety  of  similar  scenarios  for  threat 
analysis  or  system  optimization. 
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Statement  of  the  Problem 


The  objective  of  this  research  is  to  develop  and  evaluate  new 
quadrature  sets  that  produce  accurate  discrete  ordinates  solutions  to  the 
BTE  using  a  minimum  number  of  ordinates.  The  viability  of  including 
ordinates  perpendicular  to  one  or  two  cardinal  directions  is  examined.  This 
includes  deriving  and  implementing  the  appropriate  quadrature  angles  and 
weights,  and  comparing  the  results  with  those  obtained  with  standard  level- 
symmetric  LQn  quadratures  and  with  Monte  Carlo  benchmark  solutions. 

Scope 

This  research  includes  the  derivation  and  implementation  of  discrete 
ordinates  quadrature  sets  that  include  directions  perpendicular  to  one  or  two 
cardinal  directions.  Demonstration  of  the  method  including  comparison  to 
traditional  level-symmetric  quadratures  with  regard  to  computational  cost 
(execution  time)  and  accuracy  of  results  based  on  a  benchmark  calculation  is 
performed.  The  test  problems  use  three-dimensional  Cartesian  coordinates 
with  no  time  or  energy  dependence.  The  test  problems  were  run  using 
TETRAN  (13),  an  unstructured  mesh  tetrahedral  cell  code  developed  at  the 
Air  Force  Institute  of  Technology  and  the  THREEDANT  code  of  the  RSICC 
Computer  Code  Collection,  DANTSYS  3.0  (14)  from  Los  Alamos  National 
Laboratory  (LANL)  using  a  rectangular  parallelepiped  mesh.  LANL’s  MCNP 
(Monte  Carlo  Neutron  Photon)  transport  code  package  (15)  provided 


1-11 


benchmark  solutions.  Due  to  current  limitations  in  the  TETRAN  code  still  in 
development,  the  test  problems  are  defined  as  one  energy  group,  isotropic 
scatter  transport  problems.  Multiple  levels  of  spatial  mesh  refinement  are 
used.  The  method  is  tested  to  identify  any  variation  in  performance  and 
determine  an  optimal  usage.  The  scope  of  the  comparison  using  the  output  of 
the  THREEDANT  module  is  limited  due  to  the  requirement  to  modify  the 
developed  quadratures  in  order  for  the  module  to  run.  See  chapter  IV  for  a 
discussion  of  the  modifications.  No  code  changes  ore  new  modules  were 
written  to  augment  THREEDANT  to  obtain  a  more  accurate  comparison  of 
the  quadrature  sets. 

General  Approach  and  Sequence  of  Presentation 

In  chapter  II  the  integrodifferential  form  of  the  Boltzmann  transport 
equations  is  discretized  over  angle.  A  brief  discussion  of  spatial  and  energy 
discretization  is  included.  The  consequences  of  discretization  are 
enumerated.  The  method  of  generating  the  new  quadrature  sets  is  developed 
in  chapter  III.  Several  quadrature  sets  of  various  order  are  presented.  The 
method  is  implemented  using  two  test  problems  and  the  results  are 
presented  in  chapter  IV.  The  geometry  of  each  test  problem  has  been 
selected  to,  in  the  first  case,  exacerbate,  then  in  the  second,  mitigate  the 
problem  of  ray  effects.  Traditional  level  symmetric  quadratures  are  used  on 
the  same  test  problems.  Benchmark  calculations  were  performed  on  each 
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test  problem  using  a  Monte  Carlo  simulation.  The  methods  are  compared  for 
accuracy  and  computational  efficiency  and  potential  advantages  or 
disadvantages  of  the  method  identified.  Consideration  is  given  to 
smoothness,  pointwise  and  global  accuracy,  ray-effects,  and  other  systematic 
errors. 

Once  the  method  has  been  tested  and  analyzed,  recommendations  for 
use  and  for  further  research  are  given  in  the  final  chapter.  Appendices 
contain  complete  derivations  of  the  equations  used  to  generate  the 
quadratures  as  well  as  any  mathematical  routines  used  to  solve  them.  Also, 
pertinent  portions  of  input  and  output  files  of  the  test  problems  are  included 
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II.  Theory 


This  chapter  will  present  a  development  of  the  discrete  ordinates 
method,  the  criteria  for  selecting  a  quadrature,  discuss  spatial  discretization, 
and  some  the  consequences  of  applying  these  approximations. 

The  steady  state  assumptions  reduces  the  Boltzmann  transport 
equation  (Equation  1-1)  to: 

Q- V\j/(r,E,Qj  +  CT(r,E)i(/^f,E,Qj  = 

J dE'J" dQ'as(r,E  ->  E',Q •  6')i|/(r,E',6)  +  s(r,E,Q) 

In  order  to  yield  a  convenient  normalization  over  all  angles,  the  incremental 
solid  angle  is  defined  as 


do  d9sin9  _  do  dp 
2n  2  2n  2 


(H-2) 


(H-3) 


so  that 


Definitions  of  the  above  angles  are  shown  in  Figure  II- 1  (1:  11).  Equation  (II- 
1)  contains  terms  that  are  functions  of  position,  scattering  angle,  and  energy. 
The  energy  dependence  in  the  discrete  ordinates  approximation  is  most  often 
accounted  for  by  dividing  the  energy  range  of  interest  into  a  number  of 
intervals.  It  is  assumed  that  for  each  interval,  cross-sections  are  given  as 
average  values  over  the  interval  (18:  109).  The  transport  equation  is  then 


Figure  11-1:  Particle  Entering  from  Direction  Q,  Scattering  into  direction  Q' 

solved  in  each  energy  interval  as  a  mono-energetic  equation  with  particle 
contributions  scattered  from  outside  the  energy  interval  added  as  a  source 
term  and  particles  scattering  from  the  interval  of  interest  into  other  energy 
intervals  treated  as  losses.  The  resulting  equations  are  known  as  the 
multigroup  equations  (1:  61).  For  clarity,  the  remainder  of  the  derivations 
given  will  be  for  a  mono-energetic  system.  The  Boltzmann  equation  then 
becomes  the  mono-energetic  transport  equation: 
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Q-Vv|/(r,Q)  +  CT(r)v|/(r,Q)  =  JdQ'as(r,Q-Q')  v|/(r,Q')  +  s(r,Q).  (II-4) 


In  the  discrete  ordinates  method  each  term  of  the  integrodifferential 
transport  equation  is  assumed  to  be  a  separable  function  of  space  and  angle 
and  the  dependencies  are  dealt  with  separately  (3: 1-3).  Both  the  spatial  and 
angular  variables  are  required  to  be  discretized  before  the  problem  may  be 
solved  numerically.  The  angular  discretization,  being  the  focus  of  this 
research,  will  be  considered  first  and  discussed  in  some  detail  while  the 
spatial  discretization  will  be  dealt  with  later. 

Angular  Discretization 

The  majority  of  the  derivation  that  follows  is  an  extension  into  three- 
dimensions  of  the  procedure  presented  by  Lewis  and  Miller  (1).  The 
differential  scattering  cross  sections  are  expanded  in  orthogonal  Legendre 
polynomials  p^Q-Q')  where  p0  =Q  •  Q'  is  the  cosine  of  the  scattering  angle. 

The  differential  scattering  cross  section  may  be  expressed  as 


{ys(r,Q-Q')*iZ(21  +  l)asi(r)p1(Q-Q'). 


(II -5) 


The  scattering  moments  osi  are  found  from 
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asi(r)=  J  ^a8(r,n0)P!(n0),  (H-6) 

-l  z 

having  taken  advantage  of  the  orthogonality  property  of  the  Legendre 
polynomials.  The  expansion  in  Equation  (II-5)  has  been  truncated  at  L+l 
terms,  assuming  and  adequate  level  of  approximation,  taking  into 
consideration  the  degree  of  anisotropy  and  the  availability  of  cross  section 
data.  For  the  case  of  isotropic  scatter  only  the  first  term  of  Equation  (II-5)  is 
used  (L  =  0)  resulting  in 


<J,(r>6.Q')so,(f).  (II-7) 

This  derivation  will  not  assume  isotropic  scatter.  Combining  Equations  (II-4) 
and  (II- 5)  yields 


Q  •  V¥(f ,Q)  +  a(%(f ,Q)  =  I  (21  + 1) CTsi(r)J d«'Pi(6  •  Q')V(r,0')  +  s(r,Q)  (II-8) 


which  can  be  simplified  using  the  Legendre  addition  theorem  (1:  367) 


P^Q-Q') 


(II-9) 


II-4 


where  the  Ylm(Q)  are  the  spherical  harmonics  and  the  asterisks  signifies  the 

complex  conjugate.  Using  this,  the  sum  on  the  right  hand  side  of  Equation 
(II-8)  becomes 


Ld(r)  i:Yi»(n)  Jd£l’Ylm(Q')  y(r,£2').  (11-10) 


The  angular  integral  in  expression  (11-10)  is  now  just  the  coefficients 
resulting  from  the  expansion  of  the  angular  flux  in  spherical  harmonics 


L  1 

I  I  K 

1=0  m=-l 


,(?)YL(6)- 


(ii-ii) 


with 


(11-12) 


Substituting  this  into  expression  (II- 10)  gives 


I  lYim(6)asi(?)<hm(?).  (H-13) 

1=0  m=-l  V  ' 
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Replacing  the  sum  in  Equation  (II-8)  with  expression  (11-13)  yields  a  form  of 
the  transport  equation  that  is  convenient  for  discretization  in  angle 


Q  •  Vv|/(r  ,q)  +  a(r)vj;(f  ,q)  =  X  X Ybn(^)crsi  (?)|>  bn  (r)  +  s(r,6).  (11-14) 


In  the  discrete  ordinates  approximation,  Equation  (11-14)  is  enforced  only  for 
a  set  of  discrete  directions  yielding 


&n  •V\|/n(r)  +  a(r)\|/n(r)=  X  EYlm(^n)asl(r)<t)lm(?)  +  s(f>^n)>  (11-15) 


l=0m=-l 


where  i(/n(r)s  i|/(r,Qn ).  The  scalar  flux  is  approximated  by 


(j)(f)=  X  wn\)/n(r), 

n=l 


(11-16) 


and  the  flux  moments  of  Equation  (11-12)  are  approximated  by 


im 


N 

<f)=  i 

n  =  l 


wnY 


lm 


(^n)Vn(^ 


(11-17) 
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By  defining  the  right  hand  side  of  Equation  (11-15)  as  the  emission  density  or 
scattering  source  qn(?)  and  substituting  Equation  (11-17)  for  the  flux 

moments,  Equation  (11-15)  becomes 


+  n  (?)  =  0  n  (?)  (I][- 18) 


where 


qn(f)  =  s(f,nn)+  z  EY*J6nWsi(?)  s  wn-Ylm(Qn,)  ¥n'(?)  (H-19) 

V  ’  l=0m=-l  V  '  n'=l 


The  most  common  way  of  solving  Equation  (11-18)  is  the  iteration  on  the 
scattering  source  form  of  Von  Neuman’s  series  solution  (1:  80).  The  method’s 
usefulness  is  derived  from  the  fact  that  if  the  right  side  of  Equation  (11-18)  is 
known,  solving  for  i|/(r)  is  usually  straightforward.  The  iteration  is  defined 

by  modifying  Equation  (11-18)  to 


Qn-V  +  o(?)]  v)/i1+1(f)  =  q1n(r) 


(11-20) 


where  i  is  the  iteration  index.  Equation  (11-19)  likewise  becomes 
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q1n(r)  =  s(r,Qn)+  Z  ZYim(6n)<?sl(?)  Z  wn-Ylm(Qn.)  Vn'(r)*  (H-21) 

V  ’  l=Om=-l  V  '  n'=l 


The  system  of  Equations  (11-20)  and  (11-21)  must  be  solved  to  convergence  at 
each  r  of  interest. 

Before  proceeding  to  the  methodology  of  choosing  a  quadrature  set,  a 
brief  discussion  of  the  spatial  discretation  which  allows  for  the  approximation 
of  the  those  terms  which  are  now  function  of  the  spatial  variable  only  will  be 
helpful. 

Spatial  Discretization 

Discretization  of  the  spatial  variable  in  three  dimensions  can  take 
many  forms.  The  most  common  method  is  to  divide  each  of  the  three 
Cartesian  spatial  directions,  x,  y,  and  z,  into  i,  j,  and  k  intervals  respectively 
resulting  in  a  three  dimensional  grid  containing  i  x  j  x  k  rectangular 
parallelepiped  cells.  Another  method  gaining  popularity  is  to  generate  an 
unstructured  mesh  of  tetrahedra.  All  cross  sections  are  taken  to  be  piecewise 
constant  and  therefore  not  allowed  to  vary  inside  a  given  cell.  Various 
methods  exist  then  to  calculate  the  angular  flux  \\i  n  for  each  value  of  n  in 
each  cell  given  appropriate  boundary  conditions.  THREEDANT  uses 
diamond  difference  (DD)  with  negative  flux  fix  up  in  a  structured  Cartesian 
mesh,  and  TETRAN  has  the  capability  of  using  linear  characteristic  (LC), 
exponential  characteristic  (EC)  or  step  characteristic  (SC)  on  an  unstructured 
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tetrahedral  mesh.  Each  of  these  methods  and  many  more  are  discussed  and 
derived  in  detail  elsewhere  (1,2,3,8,9,16). 

Consequences  of  Discretization 

Several  problems  arise  when  using  discrete  ordinates.  Any  time 
computations  are  performed  on  a  computer,  truncation  errors  will  occur.  The 
effects  of  truncation  error  will  accumulate  with  each  computation. 

Truncation  errors  therefore  limit  the  accuracy  achievable  by  computational 
methods.  With  each  refinement  of  the  spatial  or  angular  mesh,  the 
mathematical  model  more  closely  resembles  the  continuous  analytical 
solution,  however  the  number  of  computations  also  increases.  Though 
truncation  error  associated  with  each  calculation  should  decreases  with  mesh 
refinement,  there  is  a  limit  to  the  accuracy  gained  by  continued  mesh 
refinement  as  the  number  of  calculations  gets  increasingly  large.  This  effect 
is  most  visible  in  three-dimensional  problems  where  the  number  of  cells 
increases  as  the  cube  of  the  linear  refinement.  Another,  more  systematic 
problem  arises  due  to  the  angular  discretization  called  ray  effects.  The 
phenomenon  is  most  evident  in  problems  with  localized  sources  and  small 
scattering  cross  sections  (1:195). 

When  the  scattering  cross  section  is  small,  a  substantial  percentage  of 
particles  traveling  from  a  localized  source  to  an  area  of  interest  will  be 
uncollided.  This  will  result  in  a  peaked  distribution  about  the  discrete 
ordinates.  Clearly  these  results  are  physically  unrealistic.  Ray  effects 


II-9 


manifest  as  oscillations  (1:  197),  peaks  about  a  region  where  an  ordinate  can 
be  traced  from  a  source  region  and  valleys  in  between.  As  the  order  of  the 
angular  approximation  increases  so  does  the  number  of  oscillations,  which 
also  tend  to  decrease  in  the  amount  of  deviation  from  a  smooth  curve.  If  the 
oscillations  are  uniform  about  the  correct  curve,  then  the  integral  of  the 
scalar  flux  over  the  boundary  will  still  yield  good  results  (1:  200).  When  the 
scattering  source  makes  a  large  contribution  to  the  scalar  flux,  ray  effects 
tend  to  be  mitigated.  The  scattering  source  tends  to  be  distributed  over  a 
large  area  and  is  often  nearly  isotropic.  This  gives  a  more  uniform  angular 
distribution  of  neutrons. 

Numerical  diffusion  is  a  consequence  of  the  spatial  discretization  and 
truncation  errors  due  to  spatial  differencing.  For  example,  if  a  beam  of 
neutrons  were  to  enter  the  lower  left  corner  of  a  pure  absorbing  cube  of 
material  traveling  along  the  cube  diagonal,  one  would  expect  the  attenuated 
beam  to  exit  only  at  the  upper  right  corner.  In  the  spatial  walk  of  the 
discrete  ordinates  calculation,  each  cell  is  considered  to  have  a  distribution  of 
flux  through  out  the  cell  and  flux  can  only  enter  and  exit  the  faces  of  the  cells. 
For  a  mesh  consisting  of  regular  parallelepiped  cells,  this  means  the  beam 
entering  the  cell  in  the  bottom  left  corner  will  exit  that  cell  through  the  top 
and  sides.  The  fraction  entering  the  adjacent  cells  is  determined  by  the 
incident  angle  of  the  beam.  This  will  continue  through  the  spatial  walk 
resulting  in  a  smearing  out  of  the  beam.  This  effect,  to  a  small  extent,  may 
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mitigate  ray-effects.  A  more  complete  discussion  of  numerical  diffusion  can 
be  found  in  reference  3  where  is  referred  to  as  quasi-ray  effect. 
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III.  Method 


The  primary  desire  when  selecting  a  quadrature  set  is  to  maximize 
accuracy  while  minimizing  the  number  of  ordinates.  The  primary  source  of 
computational  cost  of  obtaining  a  discrete  ordinates  solution  to  the  BTE  is  in 
the  spatial  integration.  A  complete  walk  through  the  spatial  mesh  is 
required  for  each  ordinate.  Thus  minimizing  the  number  of  ordinates 
reduces  the  computational  cost  by  minimizing  the  number  spatial  walks 
required.  Other  concerns  are  the  mitigation  of  numerical  artifacts  both 
systematic  and  nonsystematic  such  as  truncation  errors,  ray  effects,  and 
numerical  diffusion.  With  these  concerns  in  mind,  this  work  will  present  a 
method  of  generating  simultaneous  sets  of  polynomial  equations,  which 
produce  quadrature  sets  based  on  exact  integration  of  spherical  harmonics. 
Quadrature  sets  up  to  order  seven  are  presented. 

Zero  Components 

The  key  difference  in  the  quadratures  presented  here  from  those  seen 
elsewhere  is  the  addition  of  ordinates  with  zero  components  for  one  or  two  of 
the  direction  cosines.  There  have  been  several  reasons  for  not  using  zero 
components  in  the  past,  primarily  stemming  from  arguments  in  one  or  two- 
dimensional  geometry.  In  these  problems,  vertical  vector  flux  does  not 
propagate  through  the  problem.  The  infinite  path  length  resulting  from  the 
reciprocal  cosine  term  in  the  spatially  discretized  mesh  must  also  he  dealt 
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with.  Discontinuity  in  the  vector  flux  at  a  vacuum  boundary  can  lead  to 
ambiguity  in  the  meaning  of  \|/(p  =  0)  at  the  boundary  when  using  a  regular 
Cartesian  mesh.  Also,  it  can  be  shown  that  =  0)  is  not  truly  an 
independent  variable  when  discretized  in  2-D  and  therefore  contributes  little 
to  the  accuracy  of  the  solution  (19:  132). 

When  three-dimensional  quadratures  are  developed,  they  are  often  an 
extension  of  one  and  two-dimensional  cases.  The  use  of  even  order  in 
traditional  quadrature  sets,  when  applied  to  3-D  problems  has  invariably 
been  due  to  the  even  order  used  in  the  lower  dimensional  base.  The  use  of 
diamond  difference  (DD)  is  also  pervasive  in  old  computer  codes.  DD  codes 
will  not  accept  zero  components  without  special  handling  routines  being 
developed.  There  is  often  a  desire  when  developing  quadrature  sets  to 
generate  results  that  will  run  on  the  legacy  codes  without  modification.  This 
is  unfortunate,  as  modern  parallel  computers  are  capable  of  handling  far 
greater  complexities  than  the  computers  for  which  these  codes  were 
developed.  The  development  of  computer  codes  capable  of  performing 
discrete  ordinates  calculations  on  an  unstructured  mesh  eliminates  many  of 
the  problems  associated  with  special  directions.  More  accurate  characteristic 
methods  of  dealing  with  the  spatial  integration  do  not  have  the  same 
problems  dealing  with  zero  components  as  the  DD  method. 


m-2 


Quadrature  Derivation 


The  derivation  of  these  polynomial  equations  for  the  generation  of 
quadrature  sets  is  based  on  methodology  developed  by  Dr.  Kirk  Mathews  (25) 
as  an  extension  of  work  by  B.G.  Carlson  and  others  (5).  The  basic  quadrature 
sets  presented  here  are  defined  over  the  entire  unit  sphere,  but  can  be 
represented  in  the  principal  octant  or  its  edges.  The  principal  octant  is  where 
the  components  of  Q  are  all  positive.  We  require  the  quadrature  set  to  meet 
the  total  symmetry  condition.  Total  symmetry,  sometimes  referred  to  as  cubic 
symmetry,  requires  the  quadrature  set  to  remain  invariant  under  all  axis 
exchanges,  ninety-degree  rotations  about  a  cardinal  axis,  and  reflections 
across  the  x  -  y,  x  -  z,  or  y  -  z  planes.  An  axis  exchange  operation  is  the  same 
as  a  reflection  across  any  x  =  y,  x  =  -y,  x  =  z,  x  =  -z,  y  =  z,  or  y  =  -z  plane.  The 
discrete  ordinates  can  be  represented  as  points  on  the  surface  of  the  unit 
sphere.  These  points  represent  where  the  tips  of  the  unit  vectors 
corresponding  to  each  ordinate  lie  on  the  unit  sphere.  A  base  set  refers  to 
those  ordinates  in  the  original  hextant  of  the  principal  octant.  The  original 
hextant  is  a  spherical  triangle  covering  a  one  sixth  area  of  the  principal 
octant.  Figure  III-l  shows  the  unit  sphere  with  the  principle  octant  shaded 
and  an  expanded  view  of  the  principle  octant  with  the  original  hextant 
shaded.  A  complete  quadrature  set  can  be  built  by  choosing  points  in  any 
hextant  and  reflecting  the  points  by  performing  successive  axes  exchange 
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Figure  111-2:  p  -o-  %  Exchange  Operation 
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Figure  MI-3:  p  <->  r\  Exchange  Operation 


Figure  111-4:  Complete  Principal  Octant  after  %  o  r\  Exchange  Operation 
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operations  as  shown  in  Figures  III-2  through  4  (25).  The  remaining  ordinates 
are  generated  by  sequential  reflections  across  the  x-y,  x-z  and  y-z  planes. 

The  directions  Qn  constitute  a  discrete  set  of  values  of  Q  over  its  entire 
domain  (6:  2).  If  i \i  has  a  convergent  expansion  in  spherical  harmonics,  then 
from  Equations  (11-15)  and  (11-17)  we  see  that  a  necessary  condition  for 
choosing  the  ordinates  and  weights  is  that  the  spherical  harmonic 
orthogonality  condition  up  to  the  desired  order  be  satisfied,  that  is 


J"dQY}m  (q)  —  y,WnYlm  (nn  (^n  )  ~  (III-l) 

n=l 

for  all  0  <  1,1'  <  L ,  - 1  <  m  <  1 ,  and  -  T  <  m'  <  T  with  the  desired  order  of 
precision  L.  The  1  =  1’  =  m  =  m’  =  0  case  provides  the  normalization  condition 
for  the  weights, 


Zwn=l.  (HI-2) 

n 

We  also  require  thatwn  >  0  to  reduce  to  possibility  of  obtaining  physically 
unrealistic  results. 

Because  of  the  symmetry  requirements  and  the  exchangeability  of 
coordinates,  Equation  (III-l)  will  be  satisfied  if  the  moments  equations, 
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i  =  0,  1,  ...,2L  (III- 3) 


N 

J^dQ=  > 

n=l 

are  satisfied  exactly.  In  Equation  (III-3)  i  is  an  exponent  not  a  superscript. 
This  condition  is  enough  to  ensure  that  JpV^dQ  will  integrate  exactly  for 
all  0  <  i,  j,  k  <  2L  with  i  +  j  +  k  <  2L  (25).  Also  due  to  symmetry,  if  i,  j,  or  k  is 
odd,  the  integral  is  exactly  zero.  This  is  because  the  integration  in  Equation 
(III-l)  can  be  represented  as  the  integration  of  a  product  Legendre 
polynomials,  which  are  intern  polynomials  of  p  =  cos(0),  and  e1(-m  m  ^ . 
Symmetry  guarantees  the  exponential  will  ingrate  exactly  because  for  every 
cpn  there  exists  with  equal  weight  and  equal  pn  a  cpn-  that  equals  cpn  +7t. 

The  terms  in  the  summation  therefore  cancel  except  when  m  =  m’,  which 
results  in  the  exponential  reducing  to  unity.  The  odd  powered  terms  are 
exactly  zero  because  cosine  integrates  to  zero  over  the  range  of  — 1  to  1. 

Finally,  the  case  where  i  =  0  and  the  case  where  i  =  2  are  not 
independent.  For  the  i  =  0  case,  Equation  (III-3)  becomes 


J4 


N 

£2  =  1  =  w 
n=l 


n 
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and  for  the  i  =  2  case 
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(III-5) 


1  N 

J|0.2dQ  =  —  =  Xwn^n2  . 

6  n=l 

Because  of  the  symmetry  requirements,  the  p2  values  with  equal  weights 
will  always  exist  in  sets  of  three  that  sum  to  one  (this  can  readily  be  seen  by 
examining  the  most  general  case  (case  4)  in  appendix  B  with  k  =  1).  Equation 
(III-5)  therefore  becomes 


1 

3 


(III-6) 


Equation  (III-6)  and  Equation  (III-3)  are  not  linearly  independent.  We  can 
therefore  replace  the  system  of  Equations  (III-3)  with 


1 

2k +  1 


5>nHnk  . 


k=  1,  2,  L. 
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The  system  of  Equations  (III-7)  may  have  a  solution  or  solutions  providing 
the  number  of  equations,  L,  is  equal  to  the  total  number  of  free  parameters 
(degrees  of  freedom):  the  total  number  of  independent  values  for  wn  and  pn . 
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It  is  important  to  make  clear  that  the  criteria  used  for  the  selection  of 
a  quadrature  set  is  that  it  will  evaluate,  without  quadrature  error,  integrals 
of  the  sort 

for  any  N  >  1;  which  means 

fYta(n)Y,,m.(n)dn 

must  integrate  with  out  error.  This  is  why  Equation  (III-3)  requires  p1  to 
integrate  exactly  for  i  =  0  to  2L.  Symmetry  assures  getting  all  the  odd 
powers  exactly  while  the  degrees  of  freedom  in  the  weights  and  angles  are 
used  to  get  the  even  powers.  Thus  five  degrees  of  freedom  integrate  the 
coefficients  of  a  fifth  order  expansion.  P5  anisotropic  scatter  needs  this  order 
of  expansion.  Chapter  II  discusses  the  expansion  of  the  scattering  cross 
section.  With  highly  anisotropic  scattering,  a  high  order  expansion  may  be 
needed  to  get  accurate  results.  This  requires  the  discrete  ordinates  order  to 
be  at  least  as  high.  Computational  efficiency  with  highly  anisotropic 
scattering  therefore  becomes  of  even  greater  concern.  Therefore,  obtaining 
the  desired  order  of  expansion  with  fewer  ordinates  is  very  desirable. 

The  objective  is  then  to  use  the  most  reasonable  value  of  L  and 
minimize  N,  maintain  total  symmetry,  and  produce  accurate  transport 
results.  The  desired  (or  available)  precision  of  the  arithmetic,  computation 
cost,  and  the  order  of  the  Legendre  expansion  of  the  scattering  cross  section 
will  determine  the  most  reasonable  value  of  L. 
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Figure  III-4  shows  the  principle  octant  and  the  various  possible  cases 
for  points  to  be  located  and  satisfy  all  symmetry  requirements.  The  four 
cases  not  on  the  edge  of  the  octant  correspond  with  the  cases  1  through  4 
presented  by  Carlson  (5)  and  have  1,  3,  3,  and  6  points  in  the  octant  (hence  8, 
24,  24,  48  points  on  the  unit  sphere)  with  1,  2,  2,  and  3  degrees  of  freedom 
respectively.  The  three  additional  cases,  A,  B,  and  C,  are  shared  by  two  or 
more  octants.  Case  A  points  are  on  a  primary  axis.  Case  B  points  are  in  a 
zero  plane  with  the  remaining  two  components  equal.  Case  C  points  are  in  a 
zero  plane  but  the  remaining  two  components  are  not  equal.  There  are  a 
total  of  6,  12,  and  24  points  over  the  unit  sphere  with  1,  1,  and  2  degrees  of 
freedom  for  cases  A,  B,  and  C  respectively.  Table  III-2  provides  a  summary  of 
pertinent  information  regarding  each  case. 

Table  111-1:  Summary  of  Quadrature  Case  Data 


Case 

Total  Points  on 

Unit  Sphere 

Degrees  of 

Freedom 

Parameters  to 

be  Determined 

Ordinates  per 

Degree  of  Freedom 

1 

8 

1 

w 

8 

2 

24 

2 

w,p 

12 

3 

24 

2 

W,p 

12 

4 

48 

3 

w,p,r| 

16 

A 

6 

1 

w 

6 

H 

12 

1 

w 

12 

■ 

24 

2 

W,p 

12 
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Cases  A  and  1  have  the  fewest  ordinates  per  degree  of  freedom  and  are 
therefore  preferred.  Case  4  is  the  most  expensive  and  also  the  most  difficult 
to  calculate  and  is  to  be  avoided.  The  remaining  cases  are  equivalent  in  this 
regard. 

To  distinguish  these  quadrature  sets  from  those  previously  developed 
the  notation  MQn  is  adopted  for  a  quadrature  of  order  n.  To  uniquely  solve  for 
a  quadrature  set  of  order  n,  the  total  degrees  of  freedom  in  the  equation  set 
used  must  also  equal  n.  This  is  accomplished  by  selecting  a  combination  of 
cases  from  Table  III-l  such  that  the  sum  of  the  degrees  of  freedom  equals  n. 
Equation  (III-7)  is  then  used  to  determine  the  free  parameters.  Care  must  be 
taken  to  uniquely  identify  the  quadrature  being  referenced  because  there 
may  be  multiple  quadrature  sets  of  the  same  order.  For  example,  to  find  an 
MQs  quadrature  case  1  +  case  2  +  case  A  +  case  B  will  provide  the  needed 
degrees  of  freedom.  A  different  MQs  quadrature  results  from  selecting  case  4 
+  case  C.  Further  discussion  of  notation  is  presented  later. 

When  selecting  the  cases,  those  cases  with  no  free  angles  may  only  be 
selected  once;  otherwise  a  case  may  be  used  multiple  times.  Each  free 
parameter  must  be  uniquely  identified.  For  example,  for  case  3  +  case  3  + 
case  A,  the  free  parameters  would  be  {(w3 ,  P3 ,  P3 ),  (,  w3' ,  p3' ,  r|3' ),  (w^ )} 
resulting  in  an  MQ7  quadrature. 

Because  for  each  case  only  one  point  on  the  sphere  is  independent,  the 
others  resulting  from  exchange  and  reflection  operations,  we  are  able  to  write 
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a  set  of  equations  for  each  case  that  depends  only  on  the  free  parameters  of 
that  point.  This  series  of  equations  is  then  inserted  into  the  system  of 


Equations  (III -7). 


For  case  1. 


^i=,u=4l=7? 


in  the  primary  octant  with 


wj  to  be  determined  and  a  total  of  eight  points  over  the  sphere.  Therefore, 
the  contribution  of  case  1  to  Equation  (III-7)  is 


Wi 


f  1  -\2k 


\2k 


vVsJ  \s)  \fi)  u 


\2k 


\2k 


JL_ 

V3 


2k  f 
+ 


v  V3 


2k 


f  1  \2k 


V3. 


f  ^2k 


k  =  1,  2,  ...,  L. 
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The  contributions  of  the  other  cases  are  found  analogously.  A  complete 
derivation  of  each  case  is  given  in  Appendix  A  and  is  summarized  in  Table 
III-2.  Note  that  cases  2  and  3  have  the  same  form  of  contribution  equation. 

If  the  angle  solved  for  is  less  than-^=  then  this  is  a  case  2  ordinate  and  the 

V  3 

angle  is  p,  if  the  angle  is  greater  than-^=  then  a  case  3  ordinate  results  and 

V3 

the  angle  is  r| .  The  column  labeled  “Cosines”  gives  the  equations  or  values 

needed  to  find  the  initial  (p,  q,  4)  triplet.  The  rest  of  the  quadrature  is  then 
generated  by  the  operations  described  previously.  As  an  example,  for  the 
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Table  111-2:  Contribution  of  Each  Case  to  Equation  (111-7) 


MQs  quadrature  using  case  1  +  case  2  +  case  A  +  case  B,  the  following 
equations  result 
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(HI-9) 


2wa  +  8wb 


/j\k 


+  8wi 


l~\k 

^3, 


+  8w5 


2k  0i 
1^2  +  21 


f-,  2  A 

1-^2 


k  =  1,  2,  3,  4,  5. 


These  equations  were  entered  into  a  Mathematica™  (20)  notebook  and  solved 
using  built-in  functions.  A  listing  of  the  notebook  is  found  in  Appendix  B. 

The  solution  for  Equations  (III-9)  is  shown  in  Table  III-3.  Including  all 
unique  permutations  of  the  cosines  will  complete  the  principal  octant.  For 
example,  there  are  two  additional  case  A  ordinates  with  the  same  weight  as 
the  one  shown  and  having  cosines  (0,  1,  0)  and  (0,  0,  1)  respectively.  There 
are  no  additional  case  1  ordinates  as  there  are  no  more  unique  permutations 
of  the  cosines.  Figure  III-3  shows  the  distribution  of  these  points  on  the 
principal  octant.  Both  the  exact  solution  and  the  decimal  fractions  of  the 
weights  are  given.  The  weights  are  all  positive,  as  required,  which  reduces 
the  likelihood  of  physically  unrealistic  results  such  as  negative  flux.  They  are 
also  similar  in  magnitude,  providing  for  better  conditioning  of  the  problem. 
From  Table  III-l  we  see  that  there  are  fifty  total  directions  and  five  degrees 
of  freedom  for  this  quadrature.  For  comparison,  an  LQio  quadrature  also  has 
five  degrees  of  freedom  but  has  120  total  directions  over  the  unit  sphere 
requiring  2.4  times  the  computational  cost  (assuming  both  calculations 
converge  in  the  same  number  of  iterations).  The  LQ6  level  symmetric 
quadrature  has  48  directions  but  only  three  degrees  of  freedom.  Table  III-3 
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also  gives  the  base  sets  for  an  MQ3  and  an  MQ7  quadrature  set.  Table  III-4 
shows  the  number  of  discrete  ordinates  required  to  obtain  a  desired  number 
of  degrees  of  freedom  for  various  MQn  and  LQ„  quadrature  sets. 
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Table  111-3:  Some  Possible  MQn  Quadrature  Base  Sets 


weight 

1 

CO 

II 

a 

A 

1/21 

.047619047619 

1 

0 

0 

B 

9/280 

.0321428571429 

1/V2 

0 

1/V2 

1 

4/105 

.0380952380952 

1/V3 

1/V3 

1/V3 

n  =  5 

A 

4/315 

.0126984126984 

1 

0 

0 

B 

64  /  2835 
.0225749559083 

1/V2 

0 

1/V2 

1 

27  / 1280 
.0210937500000 

1/V3 

1 IS 

1/V3 

2 

14641  /  725760 
.020173335379 

3/V1T 

1  /  VTT 

l/vrr 

n  =  7 

A 

.00904818883016 

1 

0 

0 

B 

.02103246043743 

1/V2 

0 

1/V2 

C 

.00645149153857 

.954580866940172 

0 

0.29795195665031 

1 

.01827941392342 

1/V3 

1/V3 

l/V3 

2 

.01634375972737 

.875317087598172 

.34192104070871 

.34192104070871 

Table  111-4:  Comparison  of  Computational  Cost  of  MQn  versus  LQn 


Degrees 

of 

Freedom 

MQn  Cases 
Used 

LQn 

Required 

MQn:  Total 
Number  of 
Ordinates  on 
Unit  Sphere 

LQn:  Total 
Number  of 
Ordinates  on 
Unit  Sphere 

1 

A 

S2 

6 

8 

3 

A,  B,  1 

S6 

26 

48 

5 

A,  B,  1,  2 

S10 

50 

120 

7 

A,  B,  C,  1,  2 

S12/S14* 

74 

168  /  224 

9 

A,  B,  1,  2,  4 

Si6 

110 

288 

*Si2  only  has  six  degrees  ol 

:  freedom  anc 

.  S14  has  eight.  No  level  symmetric  Sn 

quadrature  has  7  degrees  of  freedom. 


For  the  MQ7  quadrature  in  Table  III-3,  the  exact  values  are  not  given, 
because  Mathematica  was  not  able  to  solve  this  set  of  equations  exactly. 
Instead  a  numerical  solutions  was  obtained  using  a  numerical  solving  built- 
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in,  NSolve.  As  the  order  of  the  quadrature  increases  the  complexity  of  the 
polynomial  system  of  equations  increases  substantially.  Only  for  the  lowest 
order  quadrature  sets  was  I  able  to  obtain  exact  solutions. 

These  quadrature  sets  were  tested  with  a  computer  algorithm  to  verify 
they  met  the  orthogonality  condition  in  Equation  (III-l)  to  double  precision. 
Many  more  combinations  of  cases  are  available  than  are  presented  here. 
Quadrature  sets  can  be  tailored  using  Table  III-2  and  Equation  3  to  best  fit 
the  geometry  of  the  problem.  Appendix  C  shows  several  additional 
quadrature  sets. 

Generating  the  sets  of  equation  to  solve  is  a  straightforward  matter, 
solving  them  is  often  a  challenge.  Some  of  the  equation  sets  do  not  have  real 
solution  that  I  have  been  able  to  find.  The  requirement  of  positive  weights 
also  eliminates  some  possibilities.  Appendix  B  contains  some  of  the  work  I 
did  in  attempting  to  get  valid  quadrature  sets.  Table  III-6  summarizes  a 
significant  number  of  possible  quadrature  set  combinations  and  the  results  of 
my  efforts.  Only  those  combinations  using  at  least  one  of  the  new  cases  are 
shown.  The  “2(3)”  notation  in  the  Cases  column  is  intended  to  show  that  the 
resulting  equations  when  selecting  a  case  2  or  a  case  3  are  mathematically 
equivalent  as  discussed  above.  An  entry  in  the  results  column  of  one  valid 
found  means  a  valid  quadrature  set  was  found  for  the  resulting  system  of 
equations,  no  valid  means  all  the  solutions  had  either  negative  weights, 
imaginary  values  or  cosines  greater  than  one.  No  solution  means  that  the 
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system  of  equations  did  not  provide  any  solutions,  good  or  bad.  Empty  spaces 
are  quadrature  combinations  I  have  not  tried  and  those  marked  with  “tried”, 
means  they  were  attempted,  but  I  was  unable  to  solve  the  equations.  The 
letters  in  the  order  column  are  a  method  to  distinguish  different  quadratures 
of  like  order,  for  example  Msb.  Valid  quadrature  sets  are  presented  in 
Appendix  C 

Table  111-3:  Case  Combination  for  Quadrature  of  Order  n 


Order 

Cases 

Results 

3 

a 

AB1 

One  valid  found 

b 

Cl 

One  valid  found 

c 

AC 

No  solution 

d 

BC 

No  solution 

~5~ 

a 

AB12(3) 

One  valid  found 

b 

AC2(3) 

One  valid  found 

c 

BC2(3) 

One  valid  found 

d 

C4  1 

One  valid  found 

e 

ABC1 

No  solution 

f 

C12(3) 

No  valid 

g 

AB4 

Tried 

h 

A23 

i 

B23 

j 

A13 

k 

B13 

~7~ 

a 

ABC  12(3) 

One  valid  found 

b 

AC23 

One  valid  found 

c 

ABC4 

Tried 

d 

AC14 

e 

AB14 

Tried 

f 

AB42(3) 

g 

AB123 

No  valid 

~ 

a 

ABC123 

Tried 

b 

ABC2(3)  4 

Tried 

c 

Etc. 
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Though  trial  and  error  can  yield  valid  quadrature  sets  it  may  be 
possible  to  gain  some  direction  when  selecting  the  cases  to  use.  Since  the 
method  is  based  on  the  exact  integration  of  spherical  harmonics,  selecting 
cases  based  on  features  of  the  spherical  harmonics  of  the  order  being  matched 
may  yield  systems  of  equations  that  are  more  readily  solved.  Figure  III-6 
(12)  shows  the  first  few  spherical  harmonics. 


Figure  111-6  :Spherical  Harmonics,  Y,m(e,<|>) 


The  three  quadratures  show  in  Table  III-3  where  tested  and  compared 
with  level  symmetric  quadratures.  These  quadratures  were  chosen  because 
they  are  the  sets  with  the  fewest  ordinates  for  the  given  order.  Also  they 
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have  ordinates  which  are  reasonable  dispersed  over  the  unit  sphere  and  the 
weights  are  all  about  the  same  order  of  magnitude.  The  MQsa  quadrature 
was  shown  in  Figure  III-5,  the  MQ3a  and  MQ7a  quadrature  sets  are  shown  in 
Figures  III- V  and  III-8. 


Figure  MI-7:  MQ3a  Quadrature  Layout 
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Figure  111-8:  MQ7a  Quadrature  Layout 
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IV.  Results 


Some  of  the  new  quadrature  sets  were  tested  on  two  problems  and 
compared  with  the  level  symmetric  LQ„  quadratures.  The  first  problem  is  a 
cube  source  inside  a  cube  shield  region  and  was  designed  to  exacerbate  ray 
effects.  The  second  problem  is  a  spherical  source  region  in  a  spherical  shield. 
This  problem  was  selected  to  reduce  ray-effects.  Due  to  the  symmetry  of  this 
problem,  the  solution  should  also  by  symmetric;  therefore  any  deviation  from 
smooth  behavior  must  be  due  to  the  method  of  calculation.  Each  test  problem 
was  run  on  the  TETRAN  code  using  the  linear  characteristic  (LC)  method  for 
handling  the  spatial  computations.  The  MQ3,  MQ5,  and  MQv  quadrature  sets 
from  Table  III-3  are  compared  with  the  standard  LQ4,  LQ10,  and  LQ16  level 
symmetric  quadrature  sets.  For  brevity,  no  letter  subscript  will  be  used  with 
the  MQn  quadratures.  Each  problem  was  also  run  on  the  THREEDANT 
module  of  the  DANTSYS  code  package  using  Diamond  Difference  with 
negative  flux  fix-up  (DZ).  Some  quadratures  with  zero  ordinates  would  cause 
the  code  to  produce  nan  (not  a  number)  values  for  some  of  the  output.  The 
case  A  ordinates  did  not  cause  problems,  but  the  case  B  and  C  would.  The 
input  file  for  THREEDANT  requires  p,  q,  and  the  weights  to  be  entered. 

The  program  calculates  the  value  of  £ .  For  the  case  B  and  C  ordinates,  when 
£  should  be  zero,  instead  a  nan  would  be  reported.  Error  handling  routines 
prevent  the  code  from  failing  but  the  output  is  not  useful.  By  reducing  the 
precision  of  the  quadrature  angles  input  into  THREEDANT,  the  code  does 
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run  and  provide  output.  Presumably  because  the  value  of  £, ,  though  very 
near  zero,  is  sufficiently  large  so  as  not  to  cause  problems. 

In  order  to  overcome  the  limitations  of  THREEDANT  with  regard  to 
the  expectation  that  the  ordinates  not  lie  on  the  edges  of  an  octant  and  to  the 
total  number  of  directions  in  a  quadrature  set,  it  was  necessary  to  make 
further  modifications.  THREEDANT  defines  a  quadrature  only  in  the 
principle  octant,  the  remaining  angles  generated  by  reflection  operations. 
This  also  requires  the  weights  that  should  lie  on  a  boundary  to  be  divided  by 
the  number  of  octants  among  which  it  is  shared.  For  example,  in  order  for 
the  normalization  of  the  weights  to  be  correct  the  weight  of  a  case  B  or  C 
ordinate  must  be  divided  by  two.  Also,  the  number  of  ordinates  in  a 
quadrature  must  match  the  number  of  ordinates  expected.  The  variable  isn 
in  the  input  file  must  be  set  to  the  value  of  n  for  the  LQn  quadrature  intended 
for  use  (there  are  provisions  for  other  quadrature  types  not  of  relevance 
here).  The  LQn  level  symmetric  quadratures  are  selected  by  default  if  no 
quadrature  information  is  given.  If  a  quadrature  set  is  entered,  the  number 
of  ordinates  supplied  for  the  principle  octant  must  match  the  number  in  a 
LQn  quadrature  for  the  value  of  n  supplied.  This  number  is  given  by 
N  =  n(n  +  2)  /  8 .  For  the  MQ3  quadrature  there  are  seven  angles  in  the 
principle  octant.  The  closest  LQn  quadrature  with  at  least  that  many 
ordinates  is  LQs  with  ten  ordinates.  Therefore,  to  run  the  MQ3  quadrature, 
the  variable  isn  must  be  set  to  eight  and  the  first  seven  ordinates  in  the 
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quadrature  array  are  filled  with  the  MQ3  values.  The  remaining  three  slots 
in  the  quadrature  array  are  fill  with  dummy  angles  with  the  weights  set  to 
zero.  This  fix  will  maintain  the  quadrature  set’s  ability  to  integrate  properly. 
The  disadvantage  of  this  fix  is  a  loss  of  computational  efficiency  because  the 
code  still  has  to  perform  computations  using  the  dummy  values  that  do  not 
contribute  to  the  solution.  Valid  comparisons  of  computational  efficiency 
between  the  MQn  and  LQn  sets  using  THREEDANT  are  therefore  limited. 

The  TETRAN  calculations  were  performed  on  the  IBM  SP  located  at 
the  major  shared  resource  center  (MSRC)  of  Wright  Patterson  AFB.  These 
problems  were  run  on  a  single  node,  consisting  of  an  RS/6000  P2SC  model 
595  processor  with  a  clock  speed  of  135MHZ  and  1Gb  of  memory  (23).  The 
THREEDANT  code  was  run  on  an  IBM  RS-6000/590  workstation  using  the 
AIX  3.2.5  operating  system  with  256  Mb  of  memory  and  operating  at  67MHz. 
The  convergence  tolerance  for  each  problem  was  10'6. 

A  benchmark  solution  was  obtained  for  each  problem  using  a  Monte 
Carlo  simulation  generated  with  MCNP  (15).  MCNP  is  a  general  three- 
dimensional  simulation  code  widely  used  and  accepted.  MCNP  results  are 
the  average  of  the  computed  quantity  over  all  histories.  Each  history  is  a 
simulation  of  one  particle’s  motion  through  the  media.  MCNP  also  provides 
the  estimated  statistical  relative  error,  R,  at  the  one  standard  deviation  level 
defined  as  the  sample  standard  deviation  divided  by  the  sample  mean. 
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Problem  Definition 


The  cube  problem  uses  the  geometry  shown  in  Figure  IV- 1.  The  mesh 
shown  is  for  clarity  and  not  used  in  calculation.  This  problem  is  a  three- 
dimensional  extension  of  the  two-dimensional  problem  presented  by  Lathrop 
(21)  and  discussed  in  references  1  and  3.  The  source  region  is  a  2x2x2  cm 


Figure  IV-1 :  Geometry  for  the  Cube  in  Cube  Problem 

cube  centered  in  a  4x4x4  cm  cube  with  the  source  normalized  to  one  over  the 
source  volume.  The  second  problem,  shown  in  Figure  IV-2,  is  a  spherical 
source  region  in  a  spherical  shield  region.  Problem  two  used  a  source  region 
with  a  radius  of  1.5  cm  normalized  to  one  over  the  volume  and  a  shield  region 
with  a  radius  of  3  cm.  Figure  IV-2  shows  the  finest  tetrahedral  mesh  used 
with  TETRAN.  The  source  and  shield  regions  have  the  same  nuclear  data 
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for  all  test  problems.  Vacuum  boundaries  are  used  outside  all  shield  regions. 


Figure  IV-2:  Geometry  for  the  Sphere  in  Sphere  Problem 


ratio  of  —  =  — .  This  results  in  a  mean  free  path  of  ^  cm.  The  MCNP 
at  3  3 

benchmark  solution  gives  the  volume  average  scalar  flux  in  each  region,  the 

average  scalar  flux  at  the  outer  surface  and  the  net  current  through  each 

surface. 

Test  Problem  One  -  Cube  Source  in  Cube  Shield 

The  first  test  problem  is  a  simple  4x4x4  cm  cube  with  vacuum 
boundaries  and  a  uniformly  embedded  isotropically-emitting  source  of 
strength  1.0  cnr3  sec1.  The  source  is  constrained  to  the  2x2x2  cm  center 
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region  of  the  cube.  The  nuclear  data  is  given  above.  The  MCNP  solution  for 
this  problem  was  run  with  one  million  histories.  The  statistical  relative 
error,  R  as  defined  above,  is  .003  for  surface  data  and  .0007  for  volume  data. 

Tetrahedral  Mesh 

Each  tetrahedral  mesh  was  generated  using  the  Pro/Mesh  module  of 
Pro/Engineer  (22).  This  problem  was  run  with  two  levels  of  mesh  refinement. 
The  meshes  are  shown  in  Figure  IV-3.  Table  IV-1  gives  the  tetrahedral  mesh 
information  for  test  problem  one.  For  clarity,  only  the  surface  cells  are  show 
for  the  fine  mesh. 


Figure  IV-3:  Tetrahedral  Meshes  Used  in  Test  Problem  One 
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Table  IV-1  :  Tetrahedral  Mesh  Data  for  Test  problem  One 


Mesh 

Total 

Tetrahedra 

Cells  in 

Source  Region 

Net  Volume  of 
Source  Region 

Average  Optical 
Thickness 

Coarse 

191 

47 

8.0 

1.203 

Fine 

1513 

162 

8.0 

0.5997 

For  a  right  regular  problem  geometry  the  tetrahedral  mesh  does  very 
well  in  matching  the  region  volumes,  even  with  a  relatively  coarse  mesh. 

Here  the  average  optical  thickness  is  defined  as  the  mean  path  length 
through  a  cell  measured  in  mean-free  paths.  The  accuracy  of  each 
quadrature  relative  to  the  benchmark  will  be  examined  first. 

Figure  IV-4  shows  contour  plots  of  the  surface  average  scalar  flux  for 
the  coarse  mesh  for  the  MQ7  and  LQie  quadrature  sets.  Figure  IV-5  shows 
the  same  data  for  the  fine  mesh.  The  plots  of  the  other  quadrature  sets  are 
very  similar  and  are  not  shown.  The  TETRAN  code  gave  the  scalar  flux  at 
the  center  of  each  cell  and  at  the  surface  of  each  cell  on  the  boundary.  To 
generate  these  plots  it  was  necessary  to  have  values  at  the  nodes.  Each  node 
value  was  approximated  as  the  average  of  the  cell  center  values  for  each 
tetrahedron  shared  by  that  node.  This  averaging  may  effect  the  variability  of 
the  actual  surface  data.  Also  from  these  figures,  the  flux  distribution’s 
dependence  on  cell  geometry  and  orientation  is  very  evident.  These  contour 
plots  show  little  dependence  on  quadrature.  Contour  plots  of  the  net  current 
out  the  cell  faces  at  the  surface  are  very  similar  to  those  for  the  scalar  flux 
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and  are  not  presented.  These  plots  indicate  that  the  new  quadratures  do  not 
cause  any  significant  loss  in  uniformity  and  that  it  is  the  spatial  mesh  which 
is  the  primary  source  of  surface  variation  in  scalar  flux  and  current.  The 
variability  of  the  data  does  decreases  slightly  with  quadrature  order.  For  low 
order  quadrature  sets  the  LQn  sets  have  less  variability  than  the  MQn  sets 
but  this  difference  is  minor.  Using  the  tetrahedral  mesh,  it  is  very  difficult  to 
get  a  qualitative  comparison  of  quadrature  sets  using  individual  data  points. 


The  integral  values  provide  a  clearer  method  of  comparison. 


Figure  IV-4:  Contour  Plot  of  the  Surface  Average  Scalar  Flux,  Cube  Problem, 


Coarse  Tetrahedral  Mesh 


Figure  IV-5:  Contour  Plot  of  Surface  Average  Scalar,  Cube  Problem,  Fine 
Tetrahedral  Mesh 

The  surface  average  scalar  flux  for  each  quadrature  and  mesh  are 
shown  in  comparison  with  the  MCNP  benchmark  calculation  in  Figures  IV-6 
and  7.  The  surface  average  is  calculated  as  the  sum  of  the  product  of  the  flux 
at  the  vacuum  interface  of  a  surface  cell  and  the  area  of  that  cell’s  face 
divided  by  the  total  surface  area, 
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The  net  current  out  of  the  cube  is  similarly  shown  in  Figures  IV-8  and  9.  The 
benchmark  is  shown  with  error  bars  indicating  the  statistical  error  R.  With 
the  exception  of  the  surface  flux  calculation  for  the  fine  tetrahedral  mesh, 
these  plot  show  there  is  little  significant  difference  in  the  accuracy  achieved 
by  the  quadratures  tested. 


Figure  IV-6:  Surface  Average  Scalar  Flux,  Cube  Problem,  Coarse  Tetrahedral 
Mesh 
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Quadrature  vs.  Surface  Average  Scalar  Flux, 
Cube  Problem,  Fine  Tetrahedral  Mesh 


Quadrature  Set 


+  MQn 
□  SQn 
- MCNP 


Figure  IV-7:  Surface  Average  Scalar  Flux,  Cube  Problem,  Fine  Tetrahedral  Mesh 


Quadrature  vs.  Net  Current  Through  Surface, 
Cube  Problem,  Coarse  Tetrahedral  Mesh 


+  MQn 
□  LQn 
- MCNP 


Quadrature  Set 


Figure  IV-8:  Net  Current  Through  the  Surface,  Cube  Problem,  Coarse 
Tetrahedral  Mesh 
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Quadrature  vs.  Net  Current  Through  the  Surface, 

Cube  Problem,  Fine  Tetrahedral  Mesh 
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Figure  IV-9:  Net  Current  Through  the  Surface,  Cube  Problem,  Fine  Tetrahedral 
Mesh 


The  absolute  value  of  the  relative  error,  e ,  is  defined  as 


benchmark  calculated  1 
4*  benchmark 


(IV- 1) 


The  absolute  value  of  the  relative  error  for  the  current  is  found  analogously 
to  the  flux.  The  flux,  current,  and  relative  error  are  summarized  in  Tables 
VI-2  and  3  for  test  problem  one.  The  value  for  s  given  for  the  MCNP  entry  is 
actually  the  statistical  error  R  and  is  shown  for  reference. 
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Table  IV-2:  Surface  Average  Scalar  Flux,  Net  Current,  and  Relative  Error,  Cube 
Problem,  Coarse  Tetrahedral  Mesh 


Quadrature 

Scalar  Flux 

8 

Net  Current 

8 

MQ3 

0.054461034 

0.05100 

0.039715347 

0.004732 

MQs 

0.054827230 

0.04462 

0.039796675 

0.002693 

MQ7 

0.054314746 

0.05355 

0.039788637 

0.002895 

LQ4 

0.054765258 

0.04570 

0.039818522 

0.002146 

LQ10 

0.055312273 

0.03617 

0.039825301 

0.001976 

LQ16 

0.055447146 

0.03382 

0.039821360 

0.002075 

MCNP 

0.05739 

0.003 

0.03990 

0.003 

Table  IV-3:  Surface  Average  Scalar  Flux,  Net  Current  and  Relative  Error,  Cube 
Problem,  Fine  Tetrahedral  Mesh 


Quadrature 

Scalar  Flux 

8 

Net  Current 

8 

MQ3 

0.054574800 

0.04902 

0.03969063 

0.005351 

mq5 

0.055480478 

0.03324 

0.03989043 

0.0003442 

MQ7 

0.054711298 

0.04664 

0.03986661 

0.0009411 

LQ4 

0.055361991 

0.03530 

0.03998596 

0.002050 

LQ10 

0.057121652 

0.004641 

0.03995183 

0.001194 

LQ16 

0.057210746 

0.003090 

0.03992792 

0.000595 

MCNP 

0.05739 

0.003 

0.03990 

0.003 

For  the  surface  values  the  LQn  quadratures  have  slightly  better  accuracy. 
The  difference  in  relative  error  between  quadrature  sets  tested  is  small  for 
the  coarse  cube.  In  the  fine  cube  the  LQn  quadrature  sets  do  better  for  scalar 
flux  but  the  MQs  and  MQ7  provide  comparable  results  for  net  current. 

Figures  IV-10  -  13  show  the  volume  average  scalar  flux  in  each  region 
and  quadrature  with  the  MCNP  solution.  The  volume  average  is  calculated 
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as  the  sum  of  the  product  of  cell  flux  and  cell  volume  divided  by  the  total 
region  volume, 


The  source  and  shield  regions  are  shown.  Tables  IV- 4  and  5  summarize  this 
data  as  well  as  present  the  relative  error. 


Quadrature  vs.  Volume  Average  Scalar  Flux, 
Shield  Region,  Cube  Problem,  Coarse  Tetrahedral  Mesh 

i 


Figure  IV-10:  Volume  Average  Scalar  Flux,  Shield  Region,  Cube  Problem  Coarse 
Tetrahedral  Mesh 
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Quadrature  vs.  Volume  Average  Scalar  Flux, 
Shield  Region,  Cube  Problem,  Fine  Tetrahedral  Mesh 


Quadrature  Set 


Figure  IV-11  Volume  Average  Scalar  Flux,  Shield  Region,  Cube  Problem,  Fine 
Tetrahedral  Mesh 


Figure  IV-12:  Volume  Average  Scalar  Flux,  Source  Region,  Cube  Problem, 


Coarse  Tetrahedral  Mesh 
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Quadrature  vs.  Volume  Average  Scalar  Flux, 
Shield  Region,  Cube  Problem,  Fine  Tetrahedral  Mesh 
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Figure  IV-13:  Volume  Average  Scalar  Flux,  Source  Region,  Fine  Cube 


Table  IV-4  :  Volume  Average  Scalar  Flux  and  Reletive  Error,  Cube  Problem, 
Coarse  Tetrahedral  Mesh 


Quadrature 

Region 

Volume  Average 
Scalar  Flux 

MCNP  Benchmark 
+/-  .07% 

Relative 

Error 

MQ3 

Source 

0.9490771 

0.93175 

0.018596 

Shield 

0.1635121 

0.16466 

0.0069711 

MQ5 

Source 

0.9269118 

0.93175 

0.0051926 

Shield 

0.1661212 

0.16466 

0.0088739 

MQ7 

Source 

0.9291961 

0.93175 

0.0027409 

Shield 

0.1658497 

0.16466 

0.0072253 

LQ4 

Source 

0.9096752 

0.93175 

0.023692 

Shield 

0.1684340 

0.16466 

0.022920 

LQ10 

Source 

0.9191761 

0.93175 

0.013495 

Shield 

0.1670303 

0.16466 

0.014395 

LQ16 

Source 

0.9209707 

0.93175 

0.011569 

Shield 

0.1668003 

0.16466 

0.012998 
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Table  IV-5:  Volume  Average  Scalar  Flux,  Fine  Mesh 


Quadrature 

Region 

Volume  Average 
Scalar  Flux 

MCNP  Benchmark 
+/-  .07% 

Relative 

Error 

MQ3 

Source 

0.9544397 

0.93175 

0.024352 

Shield 

0.1629156 

0.16466 

0.010594 

MQ5 

Source 

0.9324943 

0.93175 

0.00079882 

Shield 

0.1646805 

0.16466 

0.00012450 

MQ? 

Source 

0.9346356 

0.93175 

0.0030970 

Shield 

0.1645381 

0.16466 

0.00074031 

LQ4 

Source 

0.9160339 

0.93175 

0.016867 

Shield 

0.1663769 

0.16466 

0.010427 

LQ10 

Source 

0.9252458 

0.93175 

0.0069806 

Shield 

0.1652950 

0.16466 

0.0038564 

LQ16 

Source 

0.9269132 

0.93175 

0.0051911 

Shield 

0.1652210 

0.16466 

0.0034068 

The  MQn  quadrature  set  performed  better  than  the  LQn  quadrature 
set  for  the  volume  average  data  in  all  but  the  MQ3  case  on  the  fine  mesh. 

Computational  efficiency  was  measured  as  the  user  time  on  the  IBM 
SP  taken  to  solve  the  problem.  Table  IV-6  shows  the  user  time  for  each 
quadrature  and  each  mesh.  Recall  the  convergence  tolerance  for  all  problems 
was  lO6. 

Table  IV-6  :User  Time  in  Seconds  Taken  to  Solve  the  Cube  in  Cube  Problem 


Quadrature 

Coarse  Mesh 
(sec) 

Iterations  to 
Convergence 

Fine  Mesh 
(sec) 

Iterations  to 
Convergence 

MQ3 

9.51 

20 

96.75 

20 

MQ5 

18.43 

19 

187.55 

20 

MQ7 

26.73 

19 

274.3 

20 

LQ4 

9.58 

19 

92.77 

LQ10 

48.15 

19 

459.17 

bhhm 

LQ16 

116.75 

19 

1100.61 

20 
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For  these  quadratures  to  be  of  value,  they  must  provide  increased 
computational  efficiency  while  maintaining  comparable  error,  or  they  must 
provide  increased  accuracy  in  comparable  time.  Figures  IV- 14  through  17 
show  plots  of  user  time  versus  the  absolute  value  of  relative  error  for  the 
surface  average  scalar  flux  and  net  current  through  the  surface.  Data  points 
nearer  to  the  origin  represent  better  performance.  For  these  surface 
averaged  values  the  LQn  quadrature  sets  perform  better  on  this  problem 
except  for  the  current  through  the  surface  on  the  fine  mesh,  where  the  MQs 
and  MQ7  quadratures  perform  well.  Figures  IV- 18  and  19  show  the  volume 
average  scalar  flux  versus  relative  error.  The  MQn  quadrature  sets  have 
consistently  better  performance,  providing  less  error  for  the  computation  cost 
in  all  cases  except  for  the  MQ3  on  the  fine  mesh.  These  plots  clearly  show  the 
savings  in  computational  costs  while  maintaining  accurate  transport  results. 
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Figure  IV-14:  Computational  Efficiency,  Surface  Flux,  Cube  Problem,  Coarse 
Tetrahedral  Mesh 
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Figure  IV-15:  Computational  Efficiency,  Surface  Flux,  Cube  Problem,  Fine 
Tetrahedral  Mesh 
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User  Time  vs.  Relative  Error  Surface  Curent 
Cube  Problem,  Coarse  Tertrahedral  Mesh 
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Figure  IV-16:  Computational  Efficiency,  Surface  Current,  Cube  Problem  Coarse 
Tetrahedral  Mesh 
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Figure  IV-17:  Computational  Efficiency,  Surface  Current,  Cube  Problem  Fine 
Tetrahedral  Mesh 
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Figure  IV-18:  Computational  Efficiency,  Volume  Average  Scalar  Flux,  Cube 
Problem  Coarse  Tetrahedral  Mesh 
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Figure  IV-19:  Computational  Efficiency,  Volume  Average  Scalar  Flux,  Cube 


Problem  Fine  Tetrahedral  Mesh 


IV-21 


Parallelepiped  Mesh 

This  problem  was  run  with  one  mesh.  The  mesh  is  16x16x16  cells  in  the 
principal  octant.  Reflective  boundaries  were  used  on  three  sides.  This 
results  in  a  2x2x2  cm  cube  which  is  one  corner  of  the  test  problem  containing 
one  eighth  of  the  volume.  The  solution  to  the  remaining  seven  eighths  of  the 
problem  is  assumed  to  be  the  same  by  symmetry.  The  source  is  normalized  to 
.125  over  this  source  volume  which  is  one  eighth  of  a  complete  source  cube. 
Data  for  the  parallelepiped  mesh  is  shown  in  Table  VI-7.  The  volumes  of  the 
regions  are  easily  conserved  with  this  meshing  method.  The  number  in 
parenthesis  in  the  volume  column  is  for  an  equivalent  volume  if  reflective 
boundaries  had  not  been  used. 


Table  IV-7  :  Parallelepiped  Mesh  Data 


Total  Cells 

Cells  in  Source 

Region 

Net  Volume  of 
Source  Region 

Optical 

Thickness 

4096 

512 

1.0  (8.0) 

.09375 

Figure  IV-20  shows  a  contour  plot  of  the  cell  center  scalar  flux  for  the 
base  layer  of  cells  ( Z-  .125  cm  plane)  using  each  quadrature  set.  The  flux 
shows  a  significant  dependency  on  quadrature  for  this  mesh.  The  MQ3  case 
shows  a  very  pronounced  ray  effect.  The  higher  order  quadrature  sets  show 
better  behavior. 
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Figure  20:  Contour  Plot  of  Scalar  Flux  in  the  Z  =  ,125cm  Plane,  Cube  Problem, 
Parallelepiped  Mesh 
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Figure  IV-21:  Scalar  Flux  at  Cube  Surface 
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Figure  IV-21  (Continued):  Scalar  Flux  at  Cube  Surface 
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Figure  IV-21  shows  a  comparison  of  the  flux  at  the  top  layer  of  cells  (Z  = 
1.875  plane).  Each  successively  lower  line  represent  a  step  of  one  cell  width 
in  the  positive  y  direction  starting  at  y  =  .125.  The  lower  order  MQn 
quadratures  display  a  substantial  drop  in  flux  when  crossing  from  the  source 
to  shield  region.  This  is  still  evident  in  the  MQ7  quadrature  but  is  not  as 
pronounced.  This  is  evidently  a  ray  effect  due  to  the  ordinates  perpendicular 
to  the  axes.  Figure  IV-22  shows  the  volume  average  scalar  flux  for  each 
quadrature  compared  with  the  MCNP  benchmark.  All  quadrature  sets 
performed  well,  with  error  less  than  one  percent  except  for  LQ4.  MQ5  and 
LQ16  performed  very  well,  both  with  error  less  than  0.1  percent.  The  net 
current  through  the  surface  for  each  quadrature  and  the  MCNP  solution  are 
plotted  in  Figure  IV-23.  This  data  is  summarized  in  Table  IV-8  with  the 
relative  error  for  each  quadrature. 

Table  IV-8:  Parallelepiped  Cube  Data  Summary 


Quadrature 

Flux 

s 

Net  Current 

8 

MQ3 

0.2622893 

0.0066899 

0.059428 

0.489260 

MQ5 

0.2607804 

0.0008985 

0.059805 

0.498714 

0.2609637 

0.0016024 

0.059759 

0.497565 

LQ4 

0.2558341 

0.0180855 

0.061041 

0.529702 

LQ10 

0.2601124 

0.0016650 

0.059972 

0.502898 

LQ16 

0.2604749 

0.0002739 

0.059881 

0.500627 

Figures  IV-24  and  25  shows  the  relative  error  vs.  computation  time. 
The  lines  were  added  for  clarity  connecting  the  MQn  and  LQn  quadratures 
with  separate  lines  by  increasing  order.  When  using  the  THREEDANT  code 
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with  these  quadratures  the  convergence  time  is  not  as  direct  a  function  of  the 
number  of  angles  in  the  quadrature  set  as  it  is  when  using  TETRAN.  This  is 
primarily  due  to  the  added  dummy  angles  required.  For  Example,  both  the 
MQ3  and  MQ5  quadrature  used  isn  =  8,  requiring  the  input  of  10  ordinates. 
For  the  MQ3  case,  three  of  the  ordinates  were  dummy  values  with  zero 
weights  but  for  the  MQ5  all  ten  ordinates  are  used.  The  MQ7  quadrature 
used  isn  =  12,  requiring  21  ordinates,  five  of  which  were  dummy  values.  In 
addition,  since  THREEDANT  uses  quadrature  sets  input  by  octant, 
redundant  calculations  are  made  for  the  case  A,  B,  and  C  ordinates  that  lie 
on  octant  boundaries.  Despite  this,  the  MQn  quadratures  still  seem  to  have 
better  computational  efficiency  when  examining  integral  values.  Table  IV-9 
shows  the  computation  times  for  each  quadrature  sets. 

Table  IV-9  :  Time  by  Quadrature,  Parallelepiped  Mesh,  Cube  Probelem 


Quadrature 

Time  (sec) 

Iterations  to 
Convergence 

MQ3 

3.6 

10 

MQ5 

3.03 

8 

MQ7 

11.49 

8 

LQ4 

1.91 

9 

LQ10 

5.56 

8 

LQ16 

24.39 

8 

IV-27 


Net  Current  (n-sec^-cm'3)  c  Scalar  Flux  (n-cr 


0.2640 


_  0.2620 
CO 

E 

%  0.2600 
0) 

(/) 

I 

p  0.2580 


0.2560 


0.2540 


0.2520 


♦ 

MQn 

□ 

LQn 

_MPMP 

Quadrature  Set 


re  IV-22  :Volume  Average  Scalar  Flux,  Cube  Problem,  Parallelepiped  Mesh 


0.0650 
0.0600 
0.0550 
0.0500 
0.0450 
0.0400 
0.0350 
0.0300 

MQ3  MQ5  MQ7  LQ4  LQ10  LQ16 
Quadrature  Set 


+  MQn 
□  SQn 
- MCNP 


Figure  IV-23:  Net  Current  through  the  Surface,  Cube  Problem,  Parallelepiped 
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Figure  IV-25:  User  Time  vs.  Relative  Error  in  Net  Current  through  the  Surface, 
Cube  Problem,  Parallelepiped  Mesh 


Test  Problem  Two  -  Spherical  Source  in  Spherical  Shield 

The  second  test  problem  is  a  3  cm  radius  sphere  with  vacuum 
boundaries  and  a  uniformly  embedded  isotropically  emitting  source  of 
strength  1.0  cnr3  sec1.  The  source  is  constrained  to  a  1.5  cm  radius  sphere  in 
the  center.  The  nuclear  data  is  the  same  as  the  previous  problem.  The  same 
computer  systems  and  convergence  tolerance  that  were  used  in  problem  one 
were  used  for  problem  two. 


Tetrahedral  Mesh 

This  problem  was  run  with  three  levels  of  mesh  refinement.  The 
structure  for  each  mesh  is  shown  in  Figure  IV-26.  Table  IV- 10  gives  the 
tetrahedral  mesh  information  for  test  problem  two.  The  volumes  shown  are 
the  sums  of  the  tetrahedron  volumes  in  each  region.  This  is  compared  with 
14.14  and  98.96  cm’3  for  an  actual  spherical  volume.  For  curved  geometry  the 

Table  IV-10  :  Tetrahedral  Mesh  Data 


Mesh 

Total 

Tetrahedra 

Cells  in 

Source 

Region 

Mesh  Volume 
in  Source 
Region  (cm3) 

Mesh  Volume 
in  Shield 
Region  (cm3) 

Average 

Optical 

Thickness 

Coarse 

211 

14 

5.885 

82.24 

1.1317 

Medium 

896 

152 

11.69 

93.48 

0.84944 

Fine 

6632 

2033 

13.74 

96.46 

0.40539 

tetrahedral  mesh  does  not  conserve  volume  very  well  until  a  very  fine  mesh  is 
used.  This  is  a  fault  in  the  design  of  mesh  generators.  They  are  usually  used 
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Figure  IV-26:  Tetrahedral  Mesh  Structures  Used  in  Problem  2 
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for  finite  elements  mechanics  calculations  where  volume  is  not  an  issue. 
Figure  IV-27  shows  the  source  region  for  the  coarsest  mesh.  This  figure 
shows  how  difficult  it  is  to  mesh  curved  regions.  The  accuracy  of  each 
quadrature  will  again  be  examined  first. 


Figure  IV-28-30  shows  contour  plots  of  the  surface  average  scalar  flux 
at  the  surface  layer  of  tetrahedra  for  each  mesh  and  various  quadrature  sets. 
Due  to  the  similarity  of  the  plots,  not  all  of  the  quadratures  are  presented. 
The  same  method  of  node  averaging  was  used  to  present  this  data.  These 
contour  plots  also  show  little  dependency  on  quadrature.  To  better  show  the 
quadrature  dependence,  Figure  IV-31  shows  contour  plots  of  the  fine  sphere 
using  MQ3  and  LQ16  quadrature  sets,  looking  from  the  -x,  x,  and  0  directions 
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respectively.  The  orientation  can  be  seen  from  the  standard  arrowhead  O, 
signifying  the  positive  direction  is  out  of  the  page,  or  tail  ®,  signifying  the 

4 

positive  direction  is  into  the  page,  notation  at  the  origin  of  each  plot.  Also 
from  this  figure,  an  unusual  asymmetry  can  be  seen.  There  is  a  biasing  in 
the  positive  x  direction.  This  appears  to  be  due  to  a  problem  in  the  TETRAN 
code  that  is  currently  under  investigation.  The  degree  of  variability  in  the 
surface  average  scalar  flux  can  be  seen  in  Figure  IV-32.  This  figure  shows 
the  TETRAN  calculated  surface  flux  arbitrarily  ordered  by  magnitude  for  the 
medium  mesh.  Due  to  the  symmetry  of  the  problem  the  actual  flux  should  be 
uniform  over  the  surface.  The  shape  of  this  curve  remains  the  same  as  the 
mesh  is  made  more  or  less  refined,  however  the  magnitude  of  the  peak 
increases  to  about  .11  for  the  coarse  mesh  and  decreases  to  about  .065  for  the 
fine  mesh.  These  plots  are  not  shown.  Figures  IV-33  shows  the  net  surface 
current  in  the  same  manner  as  the  scalar  flux.  This  data  is  ordered  by  the 
magnitude  of  the  scalar  flux  to  allow  comparison  with  the  previous  figure. 
The  current  and  scalar  fluxes  have  similar  trends  but  do  vary  independently. 
Only  the  MQ7  and  LQie  data  is  presented  here,  the  other  quadratures  have 
similar  behavior. 
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Figure  IV-28:  Contour  Plot  of  Surface  Average  Scalar  Flux,  Sphere  Problem, 
Coarse  Tetrahedral  Mesh 


Figure  IV-29:  Contour  Plot  of  Surface  Average  Scalar  Flux,  Sphere  Problem, 
Medium  Tetrahedral  Mesh 
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Figure  IV-30:  Contour  Plot  of  Surface  Average  Scalar  Flux,  Sphere  Problem, 
Fine  Tetrahedral  Mesh 
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Surface  Average  Scalar  Flux 
Sphere  Problem,  Medium  Tetrahedral  Mesh 
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Figure  IV-32:  Variability  of  Scalar  Flux  at  the  Surface  of  Sphere  Problem, 
Medium  Tetrahedral  Mesh 
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Figure  IV-33:  Variability  in  Net  Current  through  the  Surface,  Sphere  Problem, 


Medium  Tetrahedral  Mesh 
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The  surface  average  scalar  flux  and  net  current  through  the  surface 
show  little  dependence  on  quadrature.  Figure  IV-34  shows  this  for  the  scalar 
flux  on  the  coarse  sphere.  The  small  differences  in  accuracy  can  be  seen  in  a 
plot  of  relative  error  versus  quadrature  as  shown  in  Figures  IV-35-37  for  the 
scalar  flux.  Plots  of  the  net  current  are  similar  and  are  not  presented. 
Regardless  of  mesh,  the  quadratures  have  comparable  performance  in 
approximating  the  surface  values.  The  relative  error  for  this  geometry  is 
much  larger  than  for  the  cube  problem  above.  This  is  primarily  attributed  to 
the  poor  job  the  mesh  does  in  matching  the  curved  surfaces.  Though  the  MQn 
sets  do  appear  to  have  less  error  for  the  surface  flux  on  the  coarser  two 
meshes,  the  difference  is  small  and,  as  can  be  seen  from  Figure  IV-34,  has 
little  significance.  This  information  is  summarized  in  Tables  VI-11-13. 

Table  IV-1 1 :  Surface  Average  Scalar  Flux  and  Net  Current,  Sphere  Problem, 
Coarse  Tetrahedral  Mesh 


Quadrature 

Scalar  Flux 

e 

Net  Current 

8 

MQ3 

0.072776076 

0.2663 

0.05724275 

0.2363 

MQ5 

0.072876460 

0.2681 

0.05722983 

0.2360 

MQ7 

0.072896331 

0.2684 

0.05722878 

0.2360 

LQ4 

0.072993465 

0.2701 

0.05720223 

0.2354 

LQ10 

0.073059529 

0.2712 

0.05723194 

0.2360 

LQ16 

0.073047543 

0.2710 

0.05723134 

0.2360 

MCNP 

0.05747 

0.003 

0.04630 

0.003 
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Figure  IV-34:  Surface  Average  Scalar  Flux,  Sphere  Problem,  Coarse  Tetrahedral 
Mesh 
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Figure  IV-35:  Relative  Error  in  Surface  Average  Scalar  Flux,  Sphere  Problem, 
Coarse  Tetrahedral  Mesh 
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Figure  IV-36:  Relative  Error  in  Surface  Average  Scalar  Flux,  Sphere  Problem, 
Medium  Tetrahedral  Mesh 
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Figure  IV-37:  Relative  Error  in  Surface  Average  Scalar  Flux,  Sphere  Problem, 
Fine  Tetrahedral  Mesh 
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Table  IV-12:  Surface  Average  Scalar  Flux  and  Net  Current,  Sphere  Problem, 
Medium  Tetrahedral  Mesh 


Quadrature 

Scalar  Flux 

s 

Net  Current 

8 

MQ3 

0.062842334 

0.0935 

0.05024186 

0.0851 

MQs 

0.062996823 

0.0961 

0.05027019 

0.0857 

MQ7 

0.062990863 

0.0960 

0.05026612 

0.0856 

LQ4 

0.063196562 

0.0996 

0.05028158 

0.0859 

LQ10 

0.063011997 

0.0964 

0.05028019 

0.0859 

LQ16 

0.063010129 

0.0964 

0.05027461 

0.0858 

MCNP 

0.05747 

0.003 

0.04630 

0.003 

Table  IV-13:  Surface  Average  Scalar  Flux  and  Net  Current,  Sphere  Problem, 
Fine  Tetrahedral  Mesh 


Quadrature 

Scalar  Flux 

8 

Net  Current 

8 

MQ3 

0.060902901 

0.0597 

0.04857312 

0.0490 

MQ5 

0.060905444 

0.0598 

0.04857166 

0.0490 

MQ7 

0.060910336 

0.0598 

0.04857333 

0.0490 

LQ4 

0.060900515 

0.0597 

0.04857391 

0.0490 

LQ10 

0.060900610 

0.0597 

0.04856642 

0.0489 

LQ16 

0.060915316 

0.0599 

0.04856830 

0.0489 

MCNP 

0.05747 

0.003 

0.04630 

0.003 

Tables  IV- 14  and  15  summarize  the  volume  average  data.  All  quadratures 
show  nearly  equal  performance  in  calculating  the  volume  average  scalar  flux 
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Table  IV- 14:  Volume  Average  Scalar  Flux  Data,  Coarse  Sphere 


Quadrature 

Region 

Volume 

Average  Scalar 
Flux 

MCNP 
Benchmark 
+/-  .2% 

s 

MQ3 

Shield 

0.259490998 

0.17601 

0.47430 

Source 

2.073454733 

1.29600 

0.59989 

MQ5 

Shield 

0.260540247 

0.17601 

0.48026 

Source 

2.059669935 

1.29600 

0.58925 

Shield 

0.260516136 

0.17601 

0.48012 

Source 

2.060075531 

1.29600 

0.58956 

LQ4 

Shield 

0.261797736 

0.17601 

0.48740 

Source 

2.043972732 

1.29600 

0.57714 

LQ10 

Shield 

0.260675794 

0.17601 

0.48103 

Source 

2.057621449 

1.29600 

0.58767 

LQ16 

Shield 

0.26067187 

0.17601 

0.48101 

Source 

2.05772263 

1.29600 

0.58775 

for  a  given  mesh.  By  examining  the  errors  in  the  source  region  shown  in 
Tables  IV- 15  and  16  an  interesting  phenomenon  can  be  detected  when 
refining  the  mesh  from  the  medium  to  the  fine  sphere.  The  error  in  the 
source  region  increases  from  about  one  percent  up  to  about  nine  percent.  The 
expected  trend  is  for  error  to  decrease  with  this  level  of  mesh  refinement. 

The  source  of  this  increase  in  error  is  suspected  to  be  the  same  unresolved 
problems  in  the  TETRAN  code  that  induced  the  biasing  mentioned  earlier. 
This  only  seems  to  become  detectable  when  running  fine  mesh  problems.  It 
is  suspected  that  a  very  small  biasing  in  the  positive  x  direction  that  only 
accumulates  significantly  for  very  fine  meshes  is  caused  by  a  tetrahedron 
splitting  algorithm  used  by  TETRAN.  The  magnitude  of  this  problem  does 
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not  appear  to  vary  with  quadrature  and  therefore  we  can  still  use  this  output 
for  relative  comparisons  between  quadrature  sets. 


Table  IV-15:  Volume  Average  Scalar  Flux  Data,  Medium  Sphere 


Quadrature 

Region 

Volume 

Average  Scalar 
Flux 

MCNP 
Benchmark 
+/-  .2% 

s 

MQ3 

Shield 

0.206680395 

0.17601 

0.17425 

Source 

1.313078744 

1.29600 

0.01318 

MQ5 

Shield 

0.20660516 

0.17601 

0.17383 

Source 

1.312631187 

1.29600 

0.01283 

MQ7 

Shield 

0.206585308 

0.17601 

0.17371 

Source 

1.312932737 

1.29600 

0.01307 

LQ4 

Shield 

0.206521942 

0.17601 

0.17335 

Source 

1.312863629 

1.29600 

0.01301 

LQ10 

Shield 

0.20670228 

0.17601 

0.17438 

Source 

1.311474207 

1.29600 

0.01194 

LQ16 

Shield 

0.206639441 

0.17601 

0.17402 

Source 

1.312184651 

1.29600 

0.01249 

Table  IV-16:  Volume  Average  Scalar  Flux  Data,  Fine  Sphere 


Quadrature 

Region 

Volume 

Average  Scalar 
Flux 

MCNP 
Benchmark 
+/-  .2% 

Absolute 
Relative  Error 

MQ3 

Shield 

0.19327010 

0.17601 

0.09806 

Source 

1.18027900 

1.29600 

0.08929 

MQ5 

Shield 

0.19331830 

0.17601 

0.09834 

Source 

1.17998870 

1.29600 

0.08951 

MQ7 

Shield 

0.19331458 

0.17601 

0.09832 

Source 

1.17996012 

1.29600 

0.08954 

LQ4 

Shield 

0.19331059 

0.17601 

0.09829 

Source 

1.17996978 

1.29600 

0.08953 

LQ10 

Shield 

0.19336766 

0.17601 

0.09862 

Source 

1.17981243 

1.29600 

0.08965 

LQ16 

Shield 

0.19334904 

0.17601 

0.09851 

- 

Source 

1.17988161 

1.29600 

0.08960 
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Computational  efficiency  was  measured  in  the  same  fashion  as  for 
problem  one.  Table  IV-17  shows  the  user  time  for  each  quadrature  and  each 
mesh. 

Table  IV-17  :User  Time  in  Seconds  Taken  to  Solve  the  Sphere  Problem 


Quadrature 

Coarse  Mesh 

Medium  Mesh 

Fine  Mesh 

MQ3 

13.46 

62.73 

531.79 

MQ5 

25.78 

121.17 

1015.44 

MQ7 

37.99 

177.87 

1494.66 

LQ4 

12.48 

58.87 

491.42 

LQ10 

62.47 

278.98 

2471.91 

LQ16 

149.05 

701.37 

5788.87 

Figures  IV-38  shows  a  plot  of  data  from  Table  IV-14  of  user  time  versus 
relative  error  in  the  volume  average  scalar  flux  for  the  coarse  mesh.  This 
curve  is  nearly  flat.  Curves  for  the  other  levels  of  mesh  refinement  and  for 
the  net  current  as  similar.  The  lack  of  significant  features  implies  the 
quadratures  have  already  converged  to  the  minimum  error  obtainable  by  the 
discrete  ordinates  method  even  for  the  lowest  order  quadrature.  This  is  not 
surprising  considering  the  simple  nature  of  the  problem. 
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Quadrature  Set 


^MQn  Source 
E  Sn  Source 
^MQn  Shield 
XLQn  Shield 


Figure  IV-38:  Relative  Error  in  Volume  Average  Scalar  Flux,  Sphere  Problem, 
Coarse  Tetrahedral  Mesh 


Parallelepiped  Mesh 

Only  one  mesh  was  used  for  this  problem.  As  in  the  cube  case, 
reflective  boundaries  were  used  on  three  sides  and  the  problem  was  run  using 
a  one  eighth  section  of  the  sphere.  The  remainder  of  the  problem  is  assumed 
to  be  the  same  by  symmetry.  Data  for  the  parallelepiped  mesh  is  shown  in 
Table  VI- 18.  The  numbers  in  parentheses  in  the  volume  columns  are  for  an 
entire  sphere  if  reflective  boundaries  had  not  been  used.  The  volumes  shown 
here  are  compared  with  98.96  and  14.14  cm3  for  an  actual  spherical  volume. 
This  method  of  mesh  generation  has  a  difficult  time  matching  curved 
surfaces. 
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Table  IV-18  :  Parallelepiped  Mesh  Data 


Total  Cells 
in  Problem 

Cells  in 
Material 

Cells  in  Source 
Region 

Net  Volume  of 
Shield  (cm3) 

Net  Volume  of 

Source  Region  (cm3) 

3375 

1464 

211 

10.02  (80.2) 

1.688  (13.5) 

The  cells  are  cubes  .1  cm  across  corresponding  to  .075  optical  thickness.  The 
source  volume  was  allowed  to  bulge  into  the  shield  region  a  small  amount  to 
allow  for  closer  volume  modeling.  Figure  IV-39  shows  contour  plots  of  the 
scalar  flux  in  the  x-y  plane  cutting  through  the  origin.  Despite  the  rough 
geometry  of  the  mesh,  the  results  are  fairly  uniform.  Figure  IV-40  shows  the 
volume  average  scalar  flux  and  Figure  41  shows  the  relative  error.  All 
quadrature  sets  did  well,  with  error  less  than  two  percent.  This  data  is 
summarized  in  Table  IV- 19 

Table  IV-19:  Sphere  Data  Summary,  Parallelepiped  Mesh 


Quadrature 

Volume  Average 
Scalar  Flux 

Relative  Error 

Net  Current 

Relative  Error 

MQ3 

0.3110153 

0.0158015 

0.047224 

0.019887 

MQ5 

0.3111058 

0.0155153 

0.047246 

0.020359 

MQ7 

0.3111027 

0.0155251 

0.047224 

0.019871 

LQ4 

0.3107919 

0.0165086 

0.047302 

0.021565 

LQ10 

0.3111415 

0.0154021 

0.047215 

0.019677 

LQ16 

0.3112716 

0.0149906 

0.047218 

0.019755 

Figure  IV-42  shows  the  net  current  by  quadrature  with  the  MCNP 
benchmark.  Again,  the  results  show  little  deviation  by  quadrature.  Figures 
IV-43  and  44  show  the  relative  error  vs.  computation  time.  Neither 
quadrature  method  seems  to  have  a  clear  advantage. 
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Table  IV-20  :  Time  by  Quadrature,  Sphere  Problem,  Parallelepiped  Mesh 


Quadrature 

Time  (sec) 

MQ3 

3.24 

MQ5 

3.3 

MQ7 

9.58 

LQ4 

2.92 

LQ10 

5.81 

LQ16 

22.97 
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0.0806736 

0.0727333 

0.0647929 

0.0568525 


10.0489121 
0.0409718 
0.0330314 
0.025091 
0.0171506 
0.00921027 


Figure  IV-39:  Comparison  of  Contour  Plots:  Scalar  Flux,  Z  =  .125  plane,  Sphere 
Problem,  Parallelepiped  mesh 
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Figure  IV-40  :  Volume  Average  Scalar  Flux,  Sphere  Problem,  Parallelepiped 
Mesh 
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Figure  IV-41:  Relative  Error,  Volume  Average  Scalar  Flux,  Sphere  Problem, 
Parallelepiped  Mesh 
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Figure  IV-42:  Net  Current,  Sphere  Problem,  Parallelepiped  Mesh 
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Figure  IV-43:  Relative  Error  in  Net  Current,  Sphere  Problem,  Parallelepiped 
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■e  IV-44:  Time  vs.  Error  in  Flux,  Sphere  Problem,  Parallelepiped  Mesh 
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Figure  IV-45:  Time  vs.  Error  in  Current,  Sphere  Problem,  Parallelepiped  Mesh 


V.  Conclusion 


The  quadratures  developed  here  show  the  potential  for  use  in  discrete 
ordinates  transport  calculations.  The  MQn  quadratures  tested  do  provide 
accuracy  comparable  to  traditional  methods,  and  often  at  a  substantial 
saving  in  time.  However,  there  is  as  yet  insufficient  data  to  make  conclusive 
statements  on  the  effectiveness  of  the  method. 

Decreased  computational  cost  for  use  on  parallelepiped  meshes  may  be 
obtained  by  development  of  new  computer  codes  to  take  advantage  of  the  one, 
or  two  dimensional  aspect  of  the  special  directions,  and  by  use  of  more 
flexible  quadrature  input  modules  to  allow  any  number  of  ordinates. 

The  MQn  quadratures  seem  to  perform  best  in  the  unstructured  mesh. 
It  is  difficult  to  determine  the  type  problem  with  regard  to  mesh  type  and 
problem  geometry  that  these  quadratures  will  work  best  with  due  to  the 
limited  amount  of  data  obtained.  The  sphere  in  sphere  problem  did  not 
adequately  differentiate  performance  between  quadrature  sets.  More  work 
is  need  on  a  greater  variety  of  problems  before  significant  conclusions  on  this 
aspect  of  performance  can  be  made. 

The  reduced  number  of  angles  in  the  new  quadratures  does  reduce  the 
computational  cost  of  the  problems  solved.  Using  the  tetrahedral  mesh,  these 
quadratures  had  comparable  accuracy  to  the  LQn  quadratures  with  regard  to 
smoothness  and  global  accuracy.  Their  performance  was  best  when  used  to 
find  integral  results.  On  the  parallelepiped  mesh,  the  MQn  quadratures 
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showed  substantial  ray  effects  though  still  did  well  when  calculating  integral 
results. 

I  recommend  using  the  MQn  quadratures  on  problems  with  little 
symmetry  on  an  unstructured  mesh.  For  structured  meshes  and  problems 
with  high  symmetry,  more  data  is  need  before  a  recommendation  can  be 
made. 


Recommendations  for  further  research 

More  quadrature  sets  need  to  be  solved  and  more  data  needs  to  be 
generated.  From  an  increased  database,  the  optimal  quadrature  set  to  use  as 
a  function  of  problem  geometry  and  mesh  type  may  be  found.  Also,  the 
potential  benefits  of  the  higher  order  of  these  quadratures  as  compared  to 
others  with  the  same  number  of  angles  needs  to  be  investigated  by  evaluating 
them  on  anisotropic  test  problems. 

Development  of  code  to  take  advantage  of  the  one-  and  two- 
dimensional  aspects  of  the  MQn  quadrature  sets  may  provide  an  interesting 
challenge  with  potentially  a  great  deal  of  gain.  Integrating  the  new  code 
modules  with  current  programs  may  allow  for  easy  transition  to  these  new 
sets  where  applicable 

I  have  been  unable  to  solve  any  quadrature  above  seventh  order. 
I  spent  many  hours  on  Mathematica,  and  attempted  to  write  a  FORTRAN 
program  using  an  expansion  method,  but  was  unsuccessful.  More  time  in 
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developing  a  method  of  solving  the  high  order  polynomials  may  yield  higher 
order  quadratures. 
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Appendix  A 


Each  case  A,  B,  C,  1,  2,  3,  and  4  included  in  a  quadrature  contribute  to  the  surrmation 

2wn^n2k  ,  k  =  l,  2,  L. 

n=l 

pua^giTec  of  the  symmety  requirements  all  of  the  weights  for  a  particular  case  are  required 
to  be  the  identical  and  *1 1  of  the  angles  for  each  case  are  functions  of  the  angles ,  u  andrj, 
in  the  base  hextant .  In  seme  of  the  cases  the  angles  have  restrictions  on  them  so  they 
are  not  free  parameters  .  Thenunberof  degrees  of  freedom  provided  by  a  case  is  the  nunber  of 
free  angles  plus  one  for  the  weight .  Each  case  below  will  have  a  discussion  of  which 
parameters  are  free . 

CaseA :  This  consist  of  six  points  over  the  unit  sphere ,  two  on  each  coordinate  axis  at  ± 

1.  All  of  the  angles  are  restricted  in  this  case  and  only  the  weight  is  free , 
yielding  one  degree  of  freedaa .  The  contribution  of  caseA.  to  the  simulation  is 

w[(l)2k+  (0)2k+  (0)2k+  (0)2k+  (0)2k+  (-l)2k] 


which  simplifies  to 


2  w 


for  all  k#  0. 

CaseB:  This  case  consist  of  12  ordiates  over  the  unit  sphere  .  Each  ordinate  has 
one  direction  cosine  equal  to  zero  and  the  other  two  equal  -  -  in  the  principal 

V  2 

octant.  The  only  free  parameter  for  this  case  then  is  the  weight  for  one  degree  of 
freedom  .  The  contribution  of  caseB  to  the  summations  is 


which  simplifies ,  fork#  0,  to 


8  w 
2k  . 


CaseC:  This  case  consist  of  24  ordinates 

over  the  unit  sphere  .  Each  ordinate  has  one  direction  cosine  equal  to  zero, 
one  free,  andtheother  is  found f ran n2  +  r?2  +  f2  =1. 

Th-ifi  leaves  one  free  angle  and  the  weight  for  two  degrees  of  freedom .  The  contribution 
to  the  summation  is 

w[4ji2k+  4  (V  1- nz  )2k+  8  (0)2k  +  4  (-ai)21c+  4  (-Vl- M2  )  ]  • 

This  simplifies  to,  fork#  0,  to 

8  w[f<2k+  ( 1  -  m2  ) k]  • 

A-l 


(  1 

\  2  k 

1  j.  A 

[  1  \ 

v  V  3 

I 

\  V  3  / 

8  w 

FT 

2k 


Casel:  This  case  consist  of  eight  ordinates  over  the  unit  sphere  .  Each  direction 
cosine  is  eqnai  to  ^  -  in  the  principal  octant .  This  leaves  only  the  weight  free 
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for  one  degree  of  freedom  .  The  contribution  is 

w[4 

This  sirrplifies  to 


Case2  :  This  consist  of  24  ordinates  over  the  unit  sphere  .  One  angle  is  free, 
the  other  two  are  equal.  This  yields  two  degrees  of  freedem, 
one  angle  and  one  weight .  The  contribution  is 


w[4/i 


2k 


8 


l-/i2 


,2k 


+  4  (-/i)2k+  8  - 


l-*2 


,2k 


] 


which  sirrplifies  to 


Case?  :  This  consist  of  24  ordinates  over  the  unit  sphere  .  One  angles  is 

free  and  the  other  two  are  equal .  This  yields  two  degrees  of  freedem  as  for  case2 
above.  There  are  two  way  to  look  at  the  contribution  for  case3  .  The  first  is  to  note 
that  all  angles  for  this  case  result  if  when  solving  for  a  case2 ,  the  angles  are  less 

than— .  This  puts  the  angle  out  of  the  base  hextant  if  this  angle  is  \x  .  But  if , 
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when  this  results ,  we  ooncider  the  angle  as  77  instead, 

this  ordinate  will  lie  in  the  base  hextant  and  is  infact  a  case3  ordinate  .  The  other 
method  is  to  proceed  as  normal  resuling  in  the  following  contribution  equation 

w[8/j21c  +  4  1-  2/l/2  j2k  +  8  (-a021c  +  4  (-V  1-  2/i2  j2  ] 


which  simplifies  to 

8w[2ji2k  +  (l-2/i2)k]  . 

I  prefer  to  use  the  former  method  and  use  a  contribution  equation  in  the  form 


The  two  equations  are  equivalent, 

as  can  be  easily  verified  by  substituting  n 


1-n2 

2 


into  the  first  result . 
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Case4  :  This  consist  of  48  directions  over  the  unit  sphere .  Both  n  and  rj  axe  free  parameters 
as  is  the  weicflit .  This  is  the  most  expensive  case  to  indued  and  the  most  difficult 
to  solve  .  To  find  the  three  angles  defining  this  ordinate  in  the  base  hextant , 
we  must  leaver  and r?  as  free  parameters  and  find  f  usings2  = 

1  -  fj2+  r)2  .  The  resulting  contribution  equation  is 

w[8^2k  +  8rj2k+  8  |Vl-^2  -  ?72)  +  8  (-m)21c+ 8  (_f7)2k+ ®  (-V 1- M2 -  H2  )  ] 

which  simplifies  to 

!6  w[jj2k+ ?}2k+  (1 -ji2  -  *72) k]  . 
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Off [General: :spelll] 
caseA  [w_,  k_]  :=  6  w  /  ;  k  ==  0 

caseA[w_,  k_]  :=  2  w  /  ;  k  !  =  0 

caseB[w_,  k_]  :  =  12  w  /  ;  k  ==  0 

caseB[w_,  k_]  :  =  8  /  ;  k  !  =  0 

caseC[w_,  ix__,  k_J  :=  24  w  /  ;  k  ==  0 

caseC[w_,  \x_,  k_]  :  =  8  w  [n2 k  +  (1  -  #i2)k)  / ;  k  !  =  0 
w 

easel  [w_,  k_]  :  =  8  — 

case2  [w_,  k_]  :  =  8  w  m2Ic  +  2 

(*  When  a  case  3  ordinate  is  desired, 
use  the  case  2  equation  an  substitute  r?  for  /i.  *) 

case4  [w_,  r)_,  k_]  :  =  16  w  (n2  k  +  (1  -  /i2  -  r?2)  +  T}2  k) 

(*  The  Filter  function  below 

searches  a  list  of  quadrature  output  from  a  Solve  function  and  *) 

(*  returns  null  values  for  elements  in  the  list  with  negative  weights  or 
angles  or  imaginary  values.  *) 

Filter  [TL_,  vars_J  :  =  Table  [If  [  (vars  /  .  TL)  [  [  j  ]  ]  =  = 

Table [Select [Re [ (vars  /  .  TL)  [  [i] ] ] ,  #  >  0&] ,  {i,  Length[TL]}] [[j]] ,  TL[[j]]] , 
{j,  Length [TL]}] 

GetQuad[Equations_,  vars_]  :=  Filter [NSolve [Equations,  vars]  ,  vars] 

1  , 

Eqns  =  Table  [Expand  [ caseA  [wA,  k]  +caseC[wC,  ixC ,  k]  ]  ==  ^  ^  ^  /  {k,  3}J 

{2  wA+  8  wC  ==  y ,  2  wA+  8  wC  -  16  wC  iiC2  +  16  wC  (uC4  ==  y ,  2  wA+  8  wC  -  24  wC  ^C2  +  24  wC  txC 

Solve  [Eqns,  {wA,  wC,  jiC}] 

{} 

(*  The  caseA  +  caseC  combination  has  no  solution  *) 

1  , 

Eqns  =  Table  [Expand  [caseB[wB,  k]  +caseC[wC,  \x  C,  k]  ]  ==  -  -  -  -  ,  {k,  3}J 

1  1 

{4  wB  +  8  wC  ==  y ,  2  wB  +  8  wC  -  16  wC  juC2  +  16  wC/jC4  ==  y,  wB  +  8  wC  -  24  wC  ^C2  +  24  wC  ^C4 


B-l 


Solve  [Eqns,  {wB,  wC,  /iC}] 

{} 

(★  The  BC  case  combination  has  no  solution  *) 

1  , 

Eqns  =  Table  [easel  [wl ,  k]  +caseC[wC,  jiC,  k]  ==  -  -  -  -  ,  {k,  3}J 

{  8^L  +  8wC==  1,  l|^  +  8wC  (MC4+  (1  -  iuC2  ) 2)  ==y,  — 2Y“  +  8  wC  (^C6  +  (1-^C2)3)  ==  y} 

Solve  [Eqns,  {wl,  wC,  /iC}] 

{{wl_>  '2§0_ '  wC^  Tlo'  13  -  V65  )  },  {wl->  jgQ-,  wC  -»  420  / 

(!3- V65  )  },  W1  ^  -sfcT'  wC^  W'  (13  +  ^5)  }' 

wC.J^,  ^^^(la  +  Vcs)}} 

(*  The  fourth  solution  above  belongs  to  the  base  set, 
the  others  coinside  with  the  ordinates 
resulting  from  symmetry  operations.  *) 

(*  Numerical 

conditioning  of  the  quadrature  will  be  best  if  the  weights  are  simmilar,  the  decimal 
values  below  show  this  to  potentialy  be  a  good  quadrature  *) 

N[%21  [  [4] ] ,  16] 

{wl  -4  0.03214285714285714,  wC  ->  0.030952380952380  95,  (uC  ->  0.9000482411921158} 

1  1 

Eqns  n  Table  [caseA[wA,  k]  +caseB[wB,  k]  +  casel[wl,  k]  ==  ,  {k,  3}J 

|8wl  +  2wa+4wB==  ~ ^  ^  ^  +  2  wA  +  2  wB  ==  — ,  ^  +  2  wA+wB  -=  y  } 

Solve  [Eqns,  {wA,  wB,  wl}] 

{{wA->  — ,  wB  -*  yQg-#  wl  280  ^ 

N[%24,  16] 

{{wA  ->  0.047  61904761904762,  wB  ->  0.0380952380952381,  wl  ->  0.03214285714285714} } 

(*  The  above  weights  are  also  similar  in  magnatude.  *) 


B-2 


1 


Eqns  =  Table  [caseA  [wA,  k]  +caseB[wB/  k]  +caseC[wC,  /jC,  k]  +  easel  [wl,  k]  =-  ^  ^  +  ^  ' 

|1^L  +  2  wA+4wB  +  8wC==  y,  +  2  wA  +  2  wB  +  8  wC  (^C4  +  (1  -  luC2)2)  =  =  y, 

■  +  2  wA  +  wB  +  8  wC  {ilCU  (1-^C2)3)  ==  y,  -yf-  +  2  wA  +  +  8  wC  (^C8  +  ( 1  -  ^C2 )  4 )  = 

■  8  wC  (fiC10  +  (1  -/JC2)5)  ==  -1-} 


3 

8  wl 
“27 

8  wl  _  _  wB 

- — —  +  2  wA  +  — ~  + 

243  4 


quadABCl  =  Solve[Eqns,  {wA,  wB,  wC,  wl,  /iC}] 


U 


(*  TheABCl  case  has  no  solutions  *) 

Eqns  =  Table [caseA[wA,  k]  +  caseB[wB,  k]  +  easel  [wl,  k]  +case2[w2,  \x2,  k]  == 


2  k+  1 


+  8  w2  +  2  wA+  4  wB  ==  y ,  -yy-  +  2  wA  +  2  wB  +  8  w2  +  y  (1  -  ^22)2)  =  = 


3 

8  wl 
27 
8  wl 
“81“ 
8  wl 
"24T 


+  2  wA  +  wB  +  8  w2  26  +  y  (1  -  }u22)  )  ==  y , 

+  2  wA  +  -y-  +  8  w2  (^28  +  y  (l-^22)4)  ==  y, 

+  2wA+^.  +  8w2  (h210+±  (1  -M22)5)  ==-b} 


quadAB12  =  Solve  [Eqns,  {wA,  wB,  wl,  w2,  /i2}] 

{{wA 

315 


4  .  64  .  27  „„  .  14641  _ 


315 

4 


2835 


1280 


/  WB  — >  o  n  o  c  /  ^1  i  oo  r\  f  ^ 


2835 


1280 


725760 

14641 

725760 


,  /^2  -» 


VTT 

3 


VTT 


■}} 


(*  The  second  solution  belongs 

to  the  base  set.  Its  numerical  value  shows  it  to  potentialy  be 
a  good  quadrature  *) 

N[quadAB12[ [2] ] ,  16] 

{wA  -»  0.0126984126984127,  wB  -»  0 . 022574 95590828924,  wl  -*  0 . 02109375, 
w2  0.02017333553791887,  tx2  -^  0.904534033733291} 


{k,  5}] 

1 

“  9  ' 


{k,  5>; 
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1 

Eqns  =  Table  [caseA  [wA,  k]  +caseC[wC,  nC,  k]  +case2[w2,  2,  k]  ==  -  —  -  -  ,  {k 


5>] 


{8  w2  +  2  wA+  8  wC  ==  y ,  2  wA  +  8  w2  (;u24  +  y  (1  -  ^22)2)  +  8  wC  (/^C4  +  (1  -  nC2)  )  ==  — , 

2wA+8w2  (aj  26+  y  (1  -  ^22)3)  +  8  wC  (;UC6  +  (1-^iC2)3)  ==  y, 

2  wA  +  8  w2  (;u28  +  y  (1  -  ^22 ) 4 )  +  8  wC  (/jC8  +  (1-juC2)4)  ==  y, 

2  wA  +  8  w2  ^ju210  +  yg-  (1  -  Ai22)  5)  +  8  wC  (,uC10  +  (1-^iC2)5)  ==  yj-} 

quadAC2  =  Solve  [Eqns,  {wA,  wC,  w2 ,  jiC,  /j2}] 


(*  Large  output  deleted  *) 

(*  I'  11  apply  the  Filter  function  to  better  see  what  I  have  *) 
Filter  [quadAC2 ,  {wA,  wC,  w2 ,  jiC,  A*2}] 


{Null,  Null,  Null,  {wA  -»  ~  V  - -f  wC-4 


2  (41  -  2  V22 ; 

2835 


44  +  13  V22 
5670 


w2  ^  11  (33-V33  (-11  +  4V22)  )  f  ^^^33-  (11  +  2  V22)  }f 


r  2  (41-2  V22  )  44  +  13  a/22 

Null,  Null,  Null,  {wA->  — i — 2835 -  '  wC  ^  - 5670 - 


w2  -4 


11  (55  -  4  V22 ) 
22680 


(33  +  ^33  (-U  +  4V22)  )  f  (11  +  2  V22)  }f  Null,  Null,  Null, 

Null,  Null,  Null,  Null,  Null} 


(*  Let1  s  see  what  the  numerical  values  are  *) 


N[%%,  16] 

{Null,  Null,  Null,  {wA-»  0.02230629169689816,  wC  -4  0 . 01851418075444525, 

w2  -4  0.017575912987  99687,  juC -»  0.507  4  563057138757,  ^2  -4  0.7858759158  676477}, 
Null,  Null,  Null,  {wA ->  0.02230629169689816,  wC  -4  0 . 01851418075444525, 
w2  -4  0.017  575912  987  99687,  ^uC  -4  0.8616774  905910132,  ^2  -4  0.7858759158  67  6477}, 
Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null} 

(*  Since  \x2  is  greater  then  this  is  a  case2  ordinate.  The  second  non 

null  solution  is  in  the  base  set.  *) 
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1 


Eqns  =  Table [caseB [wB ,  k]  +caseC[wC,  /iC,  k]  +ease2[w2,  /i2,  k]  ==  ^  ^  +  ^ 

{8  w2  +  4  wB  +  8  wC  ==  y ,  2  wB  +  8  w2  (/j24  +  y  (1  -  ^22)  )  +  8  wC  (/iC4  +  (1  -  ^C2)  )  == 
wB  +  8  w2  (/j26  +  y  (1  -  £j22)  3  j  +  8  wC  (/uC6  +  (1-^C2)3)  ==  y, 

^-  +  8w2  (^28  +  y  (1  -  ju22) 4  J  +  8  wC  (/iC8  +  (l-nC2)*)  ==  y, 

+  8  w2  (/j210  +  -i-  (1-M22)5)  +8wc  (kC10  +  (1-^C2)5)  == 


wB 


quadBC2  =  Solve  [Eqns,  {wB,  wC,  w2 ,  /iC,  A*2>] 

(*  Large  output  deleted  *) 

Filter  [quadBC2,  {wB,  wC,  w2 ,  iiC,  p2}] 

{Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null, 


{wB 

juC 


16  (-190  +  169  a/22  ) 


526995 


,  wC  -» 


169565  -  5933  a/22  11  (55  -  4  a/22  ) 


9485910 


,  w2  -4 


22680 


fjLC  -> 


/  4521  - 

U521  (2761  +  52  a/22  ) 

1 

16  ( 
\  t.yD  V  ' 

9042 

-190  +  169  a/22) 

'  wC  -» 

j  Wxj  —7 

526995  ' 

I  4521  + a} 

U521  (2761  +  52  a/22  ) 

1  9042 

,  ^2->  (11  +  2  V22)  },  Null,  Null, 


169565  -  5933  a/22  .  11  (55  -4  a/22) 

9485910  '  22680 


33 


N[%,  16] 

{Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null, 

{wB->  0.01829786661080761,  wC  -+  0.01494182037326599,  w2  -+  0.017575912987  99687, 
0.3039216446504884,  /j2  ^  0 . 7858759158676477}  ,  Null,  Null,  Null, 

{wB  -»  0.018297866610807  61,  wC  ->  0.01494182037326599,  w2  ->  0 . 01757591298799687, 
HC  -+  0.9526970315441012,  ju 2  -+  0.7858759158676477}  } 


(*  This  quadrature  also  has  promising  weights, 
the  last  solution  is  in  the  base  set.  *) 
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Eqns  =  Table [caseC [wC ,  ii£  i  k]  +  case4[w4,  /i4 ,  rj 4,  k]  ==  {k,  5)] 

|l6  w4  +  8  wC  =  =  y,  16  w4  (r?44+/j44+  (1  -  r?42  -  /j42)  2)  +  8  wC  (/;C4  +  (1-^C2)2)  ==  y, 
16  w4  (?746+^46  +  (l-r)42-/j42)3)  +  8wC  (/iC6  +  (1  ~  /^C2)  )  ==  y, 

16  w4  (r?48  +  /j48  +  (l-7742  -ju42)  )  +  8  wC  (/jC8  +  (1  -  /jC2  )  )  ==  y  r 

16  w4  (r?410  +-fj.  410  +  (1  -  r/42  -  ^42 ) 5)  +  8  wC  (fuC10  +  (1  -  /uC2)5)  == 

quadC4  =  Solve  [Eqns,  {wC,  w4 ,  fxC,  ix 4,  r/4}] 

(*  Large  output  deleted  *) 

(*  The  decimal  equivalent  after  filtering  is  shown  below  *) 

N[quadC4,  25] 
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Filter  [%,  {wC,  w4 ,  f/C,  n4 ,  rj4}] 

{Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null, 

Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null, 

Null,  {wC  -4  0.01707317073170731707317073,  -4  0 . 326499776057011143395600, 

w4  -4  0.01229674796747  967479674797,  iu4  -*  0 . 2856038324721905188708181  +  -0  .  x  10'56  I, 
r?4  -+  0. 757316019511794446172265+  -0.  x  10'57  I},  Null,  Null, 

Null,  {wC  -4  0.01707317073170731707317073,  juC  -4  0 . 32649977  6057011143395600, 
w4  ->  0.0122967  4796747  96747  9674797,  n 4  -4  0.5872843412419645862534  82  +  -0.  x 10'60  I, 
r?4  -4  0.757316019511794446172265+  -  0.  x  10'57  I},  Null,  Null, 

Null,  {wC-4  0.01707317073170731707317073,  /jC  -4  0.32649977  6057011143395600, 
w4  -+  0.01229674796747967479674797,  /u4  -+  0.587284341241964586253482+  0.  xlO'60  I, 
r)4  -+  0.2856038324721905188708181+  -  0.  x  10“5t>  1}  ,  Null,  Null, 

Null,  {wC  -4  0.01707317073170731707317073,  /iC  -4  0.326499776057011143395600, 
w4  -4  0.01229674796747  967479674797,  /J4  -4  0 . 75731601951179444 6172265 +  0 .  x  10“60  I, 

774  -4  0.2856038324721905188708181+  -0.  x  10'50  I},  Null,  Null, 

Null,  {wC  -4  0.01707317073170731707317073,  /iC  -4  0.32  64  99776057011143395600, 
w4  -4  0.01229674796747967479674797,  7j4  -4  0.2856038324721905188708181+  0.  x  10"60  I, 
774  -4  0.587284341241964586253482  +  -  0.  x  10'51  I},  Null,  Null, 

Null,  {wC  -4  0.01707317073170731707317073,  jU C -4  0.32  649977  6057011143395600, 
w4  -4  0.0122967  4796747  967479674797,  ,u4  -4  0.7573160195117944  4  6172265+  0.  xlO“55  I, 

774  -4  0.587284341241964586253482+  -  0.  x  10'51  I}, 

Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null, 

Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null, 

Null,  {wC-+  0.01707317073170731707317073,  yC  -4  0 . 945197279003024 623550787, 
w4  -4  0.01229674796747967479674797,  /u4  -4  0.2856038324721905188708181+  -0.  x  10'56  I, 
rj4  -4  0.757316019511794446172265+-0.  x  10~57  I},  Null,  Null, 

Null,  {wC->  0.01707317073170731707317073,  juC ->  0.945197279003024  623550787, 
w4  -4  0.01229674796747967479674797,  -4  0.587284341241964586253482+  -0.  xlO-60  I, 

rj4 -4  0. 757316019511794446172265+ -0.  x  10~57  I},  Null,  Null, 

Null,  {wC  -4  0.01707317073170731707317073,  nC  -4  0.94519727  9003024  623550787, 
w4  -4  0.01229674796747  967479674797,  ju4  -4  0.587284341241964586253482+  0.  x  10"60  I, 

174  -4  0.2856038324721905188708181+  -0.  x  10"50  I},  Null,  Null, 

Null,  {wC  -4  0.01707317073170731707317073,  nC  -4  0.945197279003024  623550787, 
w4  -4  0.01229674796747  9674796747  97,  ^4  -4  0.7573160195117  944461722  65+  0.  xlO'60  I, 
n4  -4  0.2856038324721905188708181+  -0.  x  10’50  I},  Null,  Null, 

Null,  {wC -»  0.01707317073170731707317073,  -4  0.94519727  9003024  623550787, 

w4  4-  0.0122967  4796747  96747  9674797,  ^4  -4  0.2856038324721905188708181+  0  .  x  10'60  I, 
r;4  -4  0.587284341241964586253482+  -  0.x  10"51  I},  Null,  Null, 

Null,  {wC->  0.017 07 317 07 317 07 31707317073,  nC  -4  0 . 945197279003024623550787, 
w4  -+  0.0122967  4796747  967479674797,  \j,4  -4  0.7573160195117  9444  61722  65+  0.  x  10’55  I, 

17 4  -4  0.5872843412419645862534  82+  -0.  x  10'51  I}} 

(*  The  filter  didi  no  eliminate  the  imaginary  portions  above  because  the 
are  artifacts  of  rounding  error  and  not  truely  imaginary  values.  Higher 
precision  math  reduces  the  imaginary  component  further.  Carfull  inspection 
shows  that  all  the  above  solutions  are  reflections  of  the  last  one, 
which  is  in  the  base  set.  *) 
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1  1 

Eqns  =  Table [easel [wl ,  k]  +  caseC[wC,  LiC,  k]  +  case2[w2,  \x2 ,  k]  ==  ^  ^  +  ”  /  {k'  ^}J 

(/j24  +  y  (1  -  /u22)2)  +  8  wC  (,uC4  +  (l-A^C2)2)  ==  1, 


1 

7  ' 


{1^L  +  8w2  +  8wC=  =  y,  l|l  +  8w2 

-^l  +  8w2  (/r26  +  y  (1-M22)3)  +8wC  (^C6  +  (1-^C2) 

4^-  +  8  w2  (^28  +  1  (1-/J22)4)  +  8  wC  (/JC8  +  (l-/iC2)4)  ==  y, 

81  \  o  > 

8  wl 

”243" 


+  8  w2  |/j210  +  -i-  (1  -V 22)5)  +  8  wC  (juC10  +  (1  -/JC2)5)  ==  -1-} 


quadlC2  =  Solve  [Eqns,  {wl ,  wC,  w2 ,  £*C,  /j2}] 
(*  Large  output  deleted  *) 


N[quadlC2 ,  40] 

Filter [%,  {wl,  wC,  w2,  juC,  Ji2}] 

{Null,  Null,  Null,  {wl-> 0.03139662 0704718662 8707375287164582 5335788+0.  xlO"52  I, 
wC -+  0.0311894  6912313230360  638905694603843725377+  0.  xlO’53  I, 
w2 -+  0.00001165730862814210336510014847547829359805+  0.  x  10“53  X, 
juC-+  0.440124599527176955730943844695059472227+  -0.  xlO'68  I, 

H2  -+  1.4828 987 31354 957 47 25183517251 8157 84 55870+  -0.  xlO’68  I},  Null, 

Null,  Null,  {wl  -+  0.03139662070471866287073752871645825335788+  0.  x 10"52  I, 
wC  -+  0.03118946912313230360638905694603843725377+  0.  xlO'53  I, 
w2  -+  0.00001165730862814210336510014847547829359805+  0.  xlO'53  I, 
luC^O.  897936710960768162219787279180959387596+  -0.  xlO'68  I, 

H2  -+  1.482898731354957472518351725181578455870+  -0.  xlO'68  I},  Null,  Null, 

Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null} 


(*  From  the  above  two  sets  of  output, 

we  see  that  there  are  no  acceptabe  1  C2  quadratures.  The  filtered  output 
shows  quadratrues  with  cosines  greater  than  one  and  from  the  unfiltered 
solutions  we  see  there  are  quadratres  with  valid  angles  but  negative  weights.  *) 

1  , 

Eqns  =  Table  [caseA  [wA,  k]  +  caseB  [wB ,  k]  + case4 [w4 ,  Ai4,  774,  k]  =  s  —  —  -  —  /  (k/  5}J 

1  2  1 
{16  w4  +  2  wA+  4  wB  ==  y,  2  wA  +  2  wB  +  16  w4  (r?44  +ju44  +  (1  -  J?42  -^42)  )  ==  y, 

2  wA  +  wB  +  16  w4  (r?46  +  /j46  +  (1  -  ??42  -  Ai42)  3)  ==  y , 

2  wA  +  -y-  +  16  w4  (r/48  +  ,u48  +  (1  -  r?42  -  /j42)  4 )  ==  y, 

2  wA  +  -y-  +  1 6  w4  (r/410  +/j410  +  (1  -  r/42  - /u42 )  5)  ==  -1-} 

quadBC2  =  Solve[Eqns,  {wA,  wB,  w4  ,  \ik ,  774}] 

(*  Large  output  deleted  *) 

N [%  ,  40] 


B-8 


Filter  [%  ,  {wA,  wB,  w4  ,  #i4,  rj 4}] 

{Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null, 
Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null,  Null} 


(*  there  are  no  valid  AB4  quadratures  *) 

(*  The  next  series  of  quadratures  are  of  order  7  *) 


Eqns  =  Table [ 

caseA[wA,  k]  +caseB[wB,  k]  +caseC[wC,  n  C,  k]  +casel[wl,  k]  +case2[w2,  \x2 ,  k] 


{ 


{k. 

7}] 

8  wl 

3 

+  8 

w2  +  2  wA 

+  4  wB 

+  8  wC  = 

1 

=  ~  3 

t 

8  wl 

9 

+  2 

wA  +  2  wB 

+  8  w2 

(/i24  + 

--M22)2 

)  +  8  wC  (ajC4  +  (1 

-MC2)2)  =  = 

8  wl 

27 

+  2 

wA  +  wB  + 

8  w2  (m26  +  1 

■  (i- 

•M22)3) 

+  8  wC  (/uC6  +  (1  - 

MC2)3)  ==l 

8  wl 
81 

+  2 

wB 

wA  +  — 

+  8  w2 

(m28  +  • 

-M22)" 

)  +  8  wC  (mC8  +  (1 

-MC2)4)  == 

8  wl 

243 

+  2 

wB 

w A  +  — r- 
4 

+  8  w2 

(m210  + 

1 

16 

(1  -M22) 

i  5)  +  8  wC  (mC10  + 

(1  -MC2)5) 

8  wl 
729 

+  2 

,  wB 
wA+  — 

+  8  w2 

(m212  + 

1 

32 

(1  -  m22) 

i6)  +  8wC  (/jC12  + 

(1  -  MC2)6) 

8  wl 

2187 

-  +  2 

>  -n  WB 

!wA+T6 

+  8  w2 

(m214  < 

1 

64 

(1-M22 

)7|  +  8  wC  (mC14  + 

(i-mc2)7) 

quadABC12  =  Solve  [Eqns,  (wA,  wB,  wC,  wl,  w2,  nC ,  fi  2}] 
(*  Large  output  deleted  *) 


1 

2  k  +  1  ' 


N[quadABC12,  16] 

Filter  [% ,  {wA,  wB,  wl,  wC,  w2,  nC,  n2}] 

{Null,  Null,  Null,  {wA ->  0 . 009048188830155413, 

wB  -»  0.0210324  6043742795,  wl  ->  0.01827941392341811,  wC  -»  0.006451491538566835, 
w2  ->  0.01634375972737  43,  ->  0.2979519566503113,  n 2  ->  0.8753170875981718}, 

Null,  Null,  Null,  {wA  ->  0.009048188830155413,  wB  ->  0 . 02103246043742795, 
wl  -»  0.01827941392341811,  wC  0.006451491538566835,  w2  ->  0.0163437597273743, 
juC->  0. 9545808669401723,  /j2  ^  0 . 8753170875981718}  ,  Null,  Null,  Null,  Null, 
Null,  Null,  Null,  Null} 


{*  The  second  of  the  above  valid  quadratures  is  in  the  base  set, 
the  other  is  a  reflection.  *) 


B-9 


1 

Eqns  =  Table  [caseA[wA,  k]  +  caseC[wC,  \x C,  k]  +  case2[w2,  \x2 ,  k]  +  case2[w3,  /j3,  k]  --  —  —  —  —  ' 
{k,  7}] 

1 8  w2  +  8  w3  +  2  wA  +  8  wC  ==  , 

2  wA  +  8  w2  (m24  +  1  (1  -m22)2)  +  8  w3  (m34  +  y  (1  -£J32)2)  +  8  wC  (mC4  +  (1  -M  C2)  )  ==  y, 

2  wA+  8  w2  (m26  +  1  (1  -  pf22 ) 3  j  +  8  w3  (m36  +  y  (1  -^(32)3)  +  8  wC  (mC6  +  (1  -^C2)  )  == 

2  wA+  8  w2  (m28  +  y  (1  -  m22)4)  +  8  w3  (m38  +  y  (1  -m32)  4J  +  8  wC  (mC8  +  (1  -/^C2)  )  ==  y , 

2  wA  +  8  w2  (m210  +  -1-  (1  -A'22)5)  +  8  w3  (m310  +  yy  (1  -  A<32)5)  +  8  wC  (/iC10  +  (1  -^C2)  )  ==  3^1 

2  wA  +  8  w2  (m212  +  -1-  (1  -  /i22)  6j  +  8  w3  (m312  +  -1-  (1  -  m32)  6)  +  8  wC  (/jC12  +  (1  -/jC2)  )  ==  yy, 

2  wA  +  8  w2  (m214  +  -1-  (1  -/i22)7)  +  8  w3  (m314  +  -1-  (1  -  m32)7)  +  8  wC  (mC14  +  (1  -/jC2)7)  ==  yy} 

quadAC23  =  Solve[Eqns,  {wA,  wC,  w2,  w3,  fiC,  m2,  fj3}] 

$Aborted 

(*  The  computer  was  not  able  to  solve  this  exactly, 

1  will  now  try  the  numerical  function  *) 

quadAC23  =  NSolve  [Eqns,  {wA,  wC,  w2 ,  w3 ,  jiC,  \x2 ,  li3}] 

$Aborted 

(*  The  computer  was  not  able  to  find  this  solution  directly  either.  X 
will  lend  some  assistance  *) 

Eqns 14  =  Table [ 

1  , 

caseA  [wA,  k]  +caseC[wC,  ft  C,  k]  +case2[w2,  \x2  ,  k]  +case2[w3,  \x2 ,  k]  ==  ,  {k,  4}J 

1 8  w2  +  8  w3  +  2  wA  +  8  wC  —  —  t 

2  wA  +  8  w2  (m24  +  y  (l-^22)2)  +  8w3  (m34  +  y  (1  -m32)2)  +  8  wC  (mC4  +  (1  -  /JC2 ) 2 )  ==  y, 

2  wA  +  8  w2  (^26  +  y  (1-M22)3)  +  8  w3  (m36  +  y  (1-M32)3)  +  8  wC  (mC6  +  (1-MC2)3)  ==  y, 

2  wA  +  8  w2  (m28  +  y  (1-M22)4)  +  8  w3  (m38  +  y  (1  -  m32  ) 4 )  +  8  wC  (mC8  +  (1-MC2)4)  ==  y} 

Solve  [Eqnsl4 ,  {wA,  wC,  w2  ,  w3}] 

{wA,  wC,  w2 ,  w3}  =  {wA,  wC,  w2,  w3}  /  .  %%[  [1]  ]  ; 

Expand  [wA]  ; 

Together [%] ; 

Cancel [ % ] ; 

PowerExpand[%]  ; 


B-10 


Simplify [%] 

(-  (1  -  2  a/C2)2  (3  +  9  a<36  -  29  //C2  +  29  a/C4  +  //32  (-1  +  54  /iC2  -  54  juC4 )  +  a/34  (- 11  -  21  iuC2  +  21  A<C4 ) )  + 
9  a/26  (-  (1  -  2  a/C2)2  + 

H34  (-39  +  210  a/C2- 210  A/C4)  +21  a/36  (1  -  5  a/C2  +  5  a/C4)  +  A<32  (19  -  105  a/C2  +  105  a/C4)  )  + 

A/22  (  (1  -  2  A/C2)2  (1  -  54  a/C2  +  54  A/C4)  + 

9  a/36  (19  -  105  a/C2 +  105  a/C4)  +  a/34  (-281  +  1716  a/C2  -  2556  a/C4  +  1680 //C6  -  840  a/C8)  + 

A/32  (109  -  653  A/C2  +  989  -  672  A/C6  +  336  juC8 )  )  - 

A/24  (  (1  -  2  a<C2)2  (-11  -  21  A/C2  +  21  a/C4)  + 

27  a/36  (13  -  70  A/C2  +  70  A/C4)  -  3  a/34  (207  -  1225  A/C2  +  1645  A/C4  -  840  A/C6  +  420  //C8)  + 

A/32  (281  -  1716  A/C2  +  2556  a/C4  -  1680  a/C6  +  840  A/C8)  ) )  / 

(630  (-1  +  a/22)  (-1  +  |z32 )  A<C2  (-1  +  //C2)  (9  a/24  A/32  (-1  +  //32)  +  (-1  +  A<32 )  (1-2  A/C2)  + 

A/22  (-9  a/34  +  (1  -  2  A/C2)2  +  4  a/32  (2  -  3  a/C2  +  3  a/C4 )  ) ) ) 

wA  =  %27 ; 

Expand [wC]  ; 

Together [%]  ; 

Cancel [%]  ; 

PowerExpand[%]  ; 

Simplify [%] 

(3  +  2  a/32  -  9  A/34  -  9  a/24  (1-18  a/32  +  21  a/34)  +  2  a/22  (1  -  54  a/32  +  81  a/34)  )  / 

(2520  A/C2  (-1  +  A/C2)  (9a/24A/32  (-1+A/32)  +  (-1  +  //32)  (l-2//C2)2  + 

A/22  (-9  a/34  +  (1-2  a/C2)2 +  4  a/32  (2  -  3  a/C2  +  3 //C4 )  ) )  ) 

wC  =  %33; 

Expand [w2] ; 

Together [%] ; 

Cancel [%3 8] ; 

PowerExpand[%] ; 


B-ll 


Simplify  [%] 

(-9M36+  (1-2ajC2)2-3^32  (5  -  24  /jC2  +  24  ^C4)  +  jU34  (23  -  84  /uC2  +  84  ^C4 ) )  / 
(630  (-1  +  ju22)  (fi22  -  ^32)  (9^24m32  (-l+^32)  +  (-l  +  *i32)  (1-2^C2)2  + 

/j22  (-9  /j34  +  (1-2  mC2)  2  +  4  ^32  (2-3  /jC2  +  3  yuC4 ) ) ) ) 

w2  =  %; 

Expand  [w3]  ; 

Together [%]  ; 

Cancel [%] ; 

PowerExpand [ % ] ; 

Simplify [%] 

(9v2e-  (1  -2iuC2)2  +a<24  (-23  +  84  ^C2  -  84  p<C4)  +  3  ju22  (5  -  24 /uC2  +  24  ,uC4 ) )  / 
(630  (/j22-/j32)  (-1  +  /J32)  (9/j24/u32  (-1  +  v22)  +  (-1+ju32)  (1-2^C2)2  + 

/j22  (-9  /j34  +  (1  -  2  ^C2)2  +  4  ,u32  (2-3  /uC2  +  3  )UC4) ) ) ) 

w3  =  % ; 

lhs5  =  caseA[wA,  5]  +  caseC[wC,  #iC,  5]  +  case2[w2,  n2 ,  5]  +case2[w3,  /i3,  5] 
Expand[%]  ; 

Together [%]  ; 

Cancel [%] ; 

PowerExpand [%]  ; 


B-12 


Simplify  [%] 


( (-23  +  25  m32  -  5  m34  +  3  M 36)  (1  -  2  /jC2)2  - 

3  /j26  ( 9  /u36  -  (1-2  mC2)2  +m34  (-23  +  84  mC2  -  84  mC4)  +  3  m32  (5  -  24  /iC2  +  24  mC4)  )  + 

iu24  (-5  (1-2  mC2 ) 2  +  m34  (65  +  636  ixC2  -  636  /jC4)  + 

3  m36  (23  -  84  mC2  +  84  mC4)  +  3  m32  (-43  -  124  mC2  +  124  /uC4)  )  + 
m22  (25  (1  -  2  nC2) 2  -  9  m36  (5  -  24  /iC2  +  24 /jC4  )  + 

m32  (149  -  112  ixC2  +  112  mC4)  +  3  ;u34  (-43  -  124  /xC2  +  124  /xC4 ) ) )  / 

(252  (9  ,u24  m32  (-1  +  m32)  +  (-1  +  m32)  (1-2/jC2)2  + 

m22  (-9  m34  +  (1-2  juC2 ) 2  +  4  m32  (2  -  3  mC2  +  3  mC4)  ) ) ) 

lhs5  =  % ; 

lhs6  =  caseA[wA,  6]  +  caseC[wC,  /iC,  6]  +  case2[w2,  Ji2 ,  6]  +case2[w3,  ix 3,  6]  ; 

Expand[%]  ; 

Together [%]  ; 

Cancel [%]; 

PowerExpand[%]  ; 

Simplify [%] 

( (1  -  2  fxC2)  2  (27  /x36  +  33  m38  + 

2  (-97-6  mC2  +  6  /UC4)  +  m32  (245  -  8  fxC2  +  8  mC4 )  -  3  m34  (37  -  12  mC2  +  12  /xC4 )  )  - 
33  ix2e  (9 /i36-  (1-2mC2)2+m34  (-23  +  84  mC2  -  84  mC4 )  +  3  m32  (5  -  24  MC2  +  24  (UC4)  )  - 
3  m26  (9 +  11  m32)  (9;u36-  ( 1  -  2  mC2 )  2  +  m34  (-23  +  84  mC2  -  84  mC4)  +  3  m32  (5  -  24  mC2  +  24  mC4 )  )  + 
M22  (  (1  -  2  mC2  )  2  (245  -  8  mC2  +  8mC4)  -  99  m38  (5  -  24  mC2  +  24  mC4  )  - 

12  m36  (31  -  151  mC2  +  151  MC4)  +  M 32  (676  +  2220  mC2  -  3948  mC4  +  3456  mC6  -  1728  mC8)  + 

54  m34  (-1  -  162  mC2  +  210  mC4  -  96  mC6  +  48  mC8)  )  +  ' 

3  m24  (-6  m36  (-7  -  6mC2  +  6mC4)  -  (1-2  fxC2)2  ( 37  -  12  mC2  +  12  mC4  )  + 

11  M38  (23  -  84  mC2  +  84  mC4)  +  18  m32  (-1  -  162  mC2  +  210  mC4  -  96  mC6  +  48  mC8)  - 

4  m34  (60  -  991  mC2  +  1243  mC4  -  504  mC6  +  252  mC8)  )  )  / 

(2520  ( 9  m24  m32  (-1+m32)  +  (-1  +  m32)  (1-2mC2)2  + 

M22  (-9  m34  +  (1-2  mC2)  2  +  4  m32  (2  -  3  mC2  +  3  mC4 )  )  )  ) 

lhs7  =  % ; 

lhs7  =  caseAfwA,  7]  +caseC[wC,  uC,  7]  +case2[w2,  m 2,  7]  +case2[w3,  M 3,  7]; 

Expand [%] ; 

Together [%]  ; 

Cancel [%]  ; 

PowerExpand[%] ; 


B-13 


Simplify [%] 

-  ( 9  ju210  (9  ju36  -  (1  -  2  iuC2)2  +1^3*  (-23  +  84  mC2  -  84  /uC4)  +  3  m32  (5  -  24  /jC2  +  24  MC4))  + 

/^28  (10  +  9/U32)  (9  /j36  -  (1-2mC2)2+m34  (-23  +  84  mC2  -  84  mC4)  +  3  m32  (5  -  24  mC2  +  24  mC4 ) )  + 
^26  (7  +  10  /j32  +  9m34) 

(9  M36  -  (1-2  MC2)2  +  m34  (-23  +  84  mC2-84mC4)  +  3  m32  (5  -  24  mC2  +  24  mC4)  )  +  (1-2mC2)2 
(47  -  7  m36  -  10  m38  -  9  m310  +  12  juC2  -  12  MC4  -  8  m32  (9  -  ju C2  +  jUC4)  +  m34  (51  -  36  mC2  +  36  MC4)  )  + 
M22  (-8  (1  -2  mC2)2  (9  -;UC2  +MC4)  +  27  m310  (5  -  24  mC2  +24  mC4)  +  3  m38  (47  -  228  MC2  +  228  mC4)  + 
M36  (95  -  464  juC2  +  464  mC4)  +  m34  (-349  +  4348  /uC2  -  6940  mC4  +  5184  mC6  -  2592  juC8)  + 

/j32  (50  -  1872  /iC2  +  3600  mC4  -  3456  iuC 6  +  1728  juC8) )  - 

M24  (-3  (1-2  iuC2 ) 2  (17-12mC2  +  12mC4)  - 

4  m36  (-5  -  24  /iC2  +  24  mC4)  +  9  /j310  (23  -  84  fu C2  +  84  mC4)  + 

M38  (95  -  192  mC2  +  192  mC4)  -  4  m34  (155  -  1376  mC2  +  2132  mC4  -  1512  mC6  +  756  mC8)  + 

m32  (349  -  4348  mC2  +  6940  mC4  -  5184  mC6  +  2592  mC8)  ) )  / 

(720  (9  m24  m32  (-1  +  m32)  +  (-1+M32)  (1-2mC2)2  + 

M22  (-9  m34  +  (1-2mC2)2  +  4m32  (2  -  3  mC2  +  3  mC4 )  )  )  ) 

lhs7  =  % ; 

1  1  1  ,  , 

NSolve[{lhs5  ==  — ,  lhs6  ==  — ,  lhs7  ==  — },  {mC,  m2,  m3}J 

$Aborted 

Solve [lhs5  ==  — ,  {/iC}] 

(*  Large  output  deleted  *) 

Length [ % ] 


(*  The  four  solutions  above  correspond  to  the  base  angle  and  its  reflections, 
I  will  continue  only  using  the  last  solution  above  *) 

jiC  =  /iC  /  .  %78  [  [4]  ] 

Expand[%]  ; 

Together [%]  ; 

PowerExpand[%]  ; 


B-14 


Simplify [%] 


(7(-l  +  23Ai32  -  55/u34  + 

33  ^36  +  33  n26  (1  -  18  ii32  +  21  ti34)  +  fi22  (23  -  448  n 32  +  1023  ii3*  -  594  n 36)  + 

11  ,u24  (-5  +  93  m32  -  159  a<34  +  63  ;u36)  - /u2  n3  V  (71  -  1807  /j32  +  8006  ^34  -  14190  ^36  + 
11187  m38  -  3267  ju310  +  1089  ju210  (-3  +  52  n32  -  18  /j34  -  204  n 3 6  +  189  ^38)  + 

33  /^28  (339  -  6287  ju32  +  9822  ju34  +  7698  /j36  -  18513  /i38  +  6237  p310)  - 
2  ju24  (-4003  +  83437  ju32  -  287774  /l(34  +  360954  ;u36  -  162063  \l3s  +  9801  /j310)  - 
22  /j26  (645  -  12706  ju32  +  32814  /j34  -  19704  /j36  -  11547  pi38  +  10098  Ii310)  + 

U22  (-1807  +  40008  ul32  -  166874  ^34  +  279532  ju36  -  207471  a<38  +  56628  u 310)  ) ) )  / 
(V2  V  (-1  +  23  ^32  -  55  ,u34  +  33  ju36  +  33  ^26  (1  -  18  n 32  +  21  ,u34)  + 

AJ22  (23  -  4  48  ju32  +  1023  /j34  -  594  ^36)  +11/J24  (-5  +  93  n32  -  159113*  +  63^36) ) ) 


MC  =  % ; 

Expand [lhs 6] ; 

Together [%]  ; 

PowerExpand[%]  ; 

Together [%]  ; 

Cancel [%]  ; 

PowerExpand[%92]  ; 

Simplify [%] 

(*  Large  ugly  equation  deleted  *) 
lhs 6  =  %; 

Expand [lhs 7] ; 

Together [%] ; 

Cancel [%] ; 

PowerExpand[%] ; 

Simplify [%] ; 
lhs7  =  % 

(*  Large  ugly  equation  deleted  *) 

NSolve[{lhs6  ==  - ,  lhs7  ==  - },  {ix2 ,  /i3}] 

L  L  13  15  J 

Out  of  memory.  Exiting. 
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quadAC23 =  FindRoot [Eqns ,  {wA,  .0130608},  {wC,  .0154866}, 

{w2,  .00673874},  {w3,  .0161761},  {/iC ,  .410254},  {^2 ,  .848421},  {fi3,  .300144}, 
AccuraeyGoal  -*  24,  WorkingPrecision  ->  34] 

{wA-*  0.01306075218457543404037831995817255,  wC ->  0 . 01548662322913343575994122509932307, 
w2  -*  0.00673874365125243712696660016440459,  w3  -4  0.01617  611174013693526966426141339587, 
0.410253515086337171560174981725413,  /. iZ  ->  0.848421498634701466179649259967932, 

M3  ->0.3001438436359286871701637504638956} 

(*  The  FindRoot  function  is  a  root  solving  function.  The  initial  guesses  used  are 
from  a  previous  effort  on  a  slower  computer.  This  is  one  possible  solution  to  the 
system  of  equations.  *) 
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Appendix  C:  Valid  Quadrature  Base  Sets 


N  =  3 


Cases 

weight 

mu 

eta 

xi 

AB 

A 

0.0476190476190476 

1 

0 

0 

B 

0.0380952380952381 

0.7071067811865475 

0 

0.7071067811865475 

1 

0.0321428571428571 

0.5773502691896258 

0.5773502691896258 

0.5773502691896258 

Cl 

c 

0.0309523809523810 

0.9000482411921158 

0 

0.4357902747044488 

1 

0.0321428571428571 

0.5773502691896258 

0.5773502691896258 

0.5773502691896258 

N  =  5 


Cases 

weight 

mu 

eta 

xi 

AB12 

A 

0.0126984126984127 

1 

0 

0 

B 

0.0225749559082892 

0.7071067811865475 

0 

0.7071067811865475 

1 

0.0210937500000000 

0.5773502691896258 

0.5773502691896258 

0.5773502691896258 

2 

0.0201733355379189 

0.9045340337332910 

0.3015113445777637 

0.3015113445777637 

AC2 

A 

0.0223062916968982 

1 

0 

0 

C 

0.0185141807544452 

0.8616774905910132 

0 

0.5074563057138757 

2 

0.0175759129879969 

0.7858759158676477 

0.4372636760921183 

0.4372636760921183 

B 

0.0182978666108076 

0.7071067811865475 

0 

0.7071067811865475 

C 

0.0149418203732660 

0.9526970315441012 

0 

0.3039216446504885 

2 

0.0175759129879969 

0.7858759158676477 

0.4372636760921183 

0.4372636760921183 

C4 

C 

0.0170731707317073 

0.9451972790030246 

0 

0.3264997760570111 

4 

.01229674796747967 

0.7573160195117944 

0.5872843412419646 

0.2856038324721905 

N  =  7 


Cases 

weight 

mu 

eta 

xi 

ABC12 

A 

0.0090481888301554 

1 

0 

0 

C 

0.0064514915385668 

0.9545808669401723 

0 

0.2979519566503114 

B 

0.0210324604374280 

0.7071067811865475 

0 

0.7071067811865475 

1 

0.0182794139234181 

0.5773502691896258 

0.5773502691896258 

0.5773502691896258 

2 

0.0163437597273743 

0.8753170875981718 

0.7854361833270270 

0.7854361833270270 

AC23 

A 

0.0130607521845754 

1 

0 

0 

C 

0.0154866232291334 

0.911971520037388 

0 

0.4102535150863372 

2 

0.0067387436512524 

0.8484214986347015 

0.374286628571238 

0.374286628571238 

3 

0.0161761117401369 

0.3001438436359287 

0.674504882535127 

0.674504882535127 
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wgt=0. 0022 620472075388518  0.0032257457692834124  0.0163437597273742918 

0.0032257457692834124  0.010516230218714844  0.01827941392341814 

0.010516230218714844  0.0032257457692834124  0.0163437597273742918 

0.0163437597273742918  0.0032257457692834124  0.0022620472075388518 


42.  0.0032257457692834124  0.010516230218714844  0.0032257457692834124 

43.  0.0022620472075388518  00000 

44.  mu=0. 0000000000000000000  0.297951956650311267  0.34192104071 
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